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Abstract 

The  optimization  of  the  three  gains  of  a  third-order 
baro-inertial  vertical  channel  has  been  formulated  as  a 
stochastic  optimal  control  problem,  with  the  objective 
of  minimizing  the  mean  squared  altitude  error  due  to  the 
noise  induced  altitude  error  and  a  disturbance  of  known 
magnitude . 

For  a  vehicle  carrying  out  a  TERCOM-update  immediately 
following  a  vertical  descent,  and  being  subjected  to  a 
disturbance  input  to  the  vertical  channel,  optimum  gains 
are  presented  and  the  performance  is  analyzed  through  a 
simulated  flight  in  a  Monte  Carlo  analysis.  Performance 
comparisons  between  the  optimized  gains  and  the  classical 
gains  are  also  presented.  The  results  show  a  significant 
performance  improvement  over  the  classical  gains  for  a 
vehicle  carrying  out  the  TERCOM-update. 
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INVESTIGATION  OF  A  THIRD -ORDER 
BARO-DAMPED  VERTICAL  CHANNEL  OF  INS 

I .  Introduction 

Background 

The  development  of  highly  accurate,  self-contained 
inertial  navigation  systems  (INS)  has  been  one  of  the 
major  engineering  accomplishments  of  the  past  fifty  years. 
In  the  simplest  terms,  an  inertial  navigation  system  is 
one  which  uses  Newton's  law  of  motions  and  a  set  of 
initial  conditions  to  determine  continuously  the  velocity, 
position  and  attitude  of  the  vehicle  in  which  it  is  con¬ 
tained.  The  first  aircraft  navigation  systems  were  pri¬ 
marily  two-channel  systems  that  provided  horizontal 
navigation  data  (Refs  1,2,3).  Inertial  navigators  using 
three  channels  were  introduced  with  the  advent  of  the 
missile  and  space  era.  In  addition,  the  value  of 
inertially  derived  altitude  and  vertical  velocity  was 
recognized  in  aircraft  and  missile  applications  involving 
low-level  flights  and  precision  weapon  delivery  (Ref  4) . 

Recently,  the  vertical  channel  performance  has 
become  important  for  a  different  reason.  Long  flights 
require  some  navigation  update  to  counter  the  long  term 
INS  drifts.  For  cruise  missiles,  a  position  update  has 
been  developed  based  on  pattern  recognition  of  the  terrain 


altitude  profile.  To  measure  this  altitude  profile,  a 
radar  altimeter  measures  terrain  clearance  while  a  constant 
altitude  flight  path  is  maintained.  The  constant  altitude 
flight  path  depends  on  the  indicated  altitude  from  the 
INS.  Clearly,  INS  errors  directly  corrupt  the  TERCOM 
(Terrain  Contour  Mapping)  data.  For  this  reason,  accurate 
altitude  tracking  by  the  INS  is  of  critical  importance 
during  the  data  taking  period  of  TERCOM. 

If  one  analyzes  the  error  behavior  of  a  local-level 
inertial  navigation  system,  one  finds  that,  given  the 
Schuler,  the  Foucault  and  the  24-hour  oscillations,  the 
horizontal  axes  (east,  north)  display  a  stable  navigation 
error  behavior,  while  the  vertical  channel  is  unstable 
(Ref  5);  that  is  to  say,  the  vertical  velocity  and  posi¬ 
tion  errors  increase  exponentially  with  the  passage  of 
time  (Appendix  A) .  This  instability  is  due  to  the  calcu¬ 
lation  of  the  gravity  correction;  an  accelerometer  measures 
all  the  accelerations  to  which  the  vehicle  is  subjected 
with  the  exception  of  acceleration  due  to  gravity.  When 
INS  acceleration  is  estimated  by  adding  measured  specific 
force  to  gravity  computed  from  a  gravity  model,  an  error 
feedback  is  established  due  to  evaluating  the  gravity  model 
with  an  imperfect  position  estimate.  In  the  vertical 
channel,  this  feedback  is  positive;  that  is,  a  positive 
vertical  position  error  creates  a  positive  vertical  accel¬ 
eration  error.  The  essentially  unstable  nature  of  such  a 


vertical  channel  mechanization  results  in  errors  which 
grow  exponentially  with  an  approximate  ten  minute  time 
constant  (Ref  5).  Thus,  for  a  typical  navigation  flight, 
the  vertical  channel  needs  to  be  stabilized  by  some  exter¬ 
nal  altitude  reference,  usually  barometric  altimeter  data. 
Unfortunately,  as  normally  implemented,  this  method  has 
one  small  drawback.  The  time  constant  associated  with  a 
barometric  altimeter  is  very  large,  which  means  that 
through  prolonged  descent  or  turns,  the  vertical  channel 
inherits  an  error  which  can  persist  for  as  long  as  two 
minutes.  This  error  can  degrade  weapon  delivery,  especially 
upon  reattack. 

A  classical  approach  to  improving  and  stabilizing 
the  vertical  channel  is  to  introduce  external  altitude 
information  from,  for  example,  a  barometric  altimeter. 

The  baro-damped  vertical  mechanization  has  evolved  to  a 
"third-order"  mechanization  which  feeds  back  two  terms  to 
the  vertical  acceleration  calculation,  and  one  term  to 
the  vertical  velocity  calculation.  The  basic  difference 
between  baro-altimeter  and  the  INS  altiutde  is  fed  back 
to  the  velocity  calculation  with  a  gain  of  .  This 
difference  is  also  fed  back  to  the  acceleration  calcula¬ 
tion  with  a  proportional  gain  of  I<2  and  an  integral  gain 
of  .  In  this  manner,  stable  vertical  channel  operation 
has  been  developed  which  has  proven  acceptable  for  many 
applications. 
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Unfortunately,  the  classical  gains  (K^f  K 2  and  K.^) 
result  in  a  sluggish  response  to  low  frequency  baro- 
induced  altitude  disturbances  which  are  encountered  during 
prolonged  descents  as  might  precede  a  TERCOM-update . 

Recent  research  in  vertical  velocity  improvements  (Ref  4) 
suggests  that  vertical  position  estimates  might  be  simi¬ 
larly  improved  based  on  optimizing  these  gains. 

Problem 

This  thesis  addresses  the  task  of  optimizing  the 
vertical  position  estimates  of  a  baro-damped  INS.  Speci¬ 
fically,  the  third-order  baro-damped  system  is  treated  to 
optimize  the  transient  vertical  performance  by  selecting 
proper  gains  (K^,  and  for  a  third-order  mechaniza¬ 
tion)  during  a  TERCOM-type  (Terrain  Contour  Mapping) 
update  following  specific  disturbance  profile  to  the  ver¬ 
tical  loop  caused  by  vehicle  maneuvers  (horizontal  or 
vertical  turns) . 

Objectives 

The  objectives  are  to  calculate  optimal  gains  for 
the  stated  problem,  investigate  sensitivities  of  per¬ 
formance  to  these  gains,  and  to  validate  the  optimal 
gains  in  a  Monte  Carlo  study. 


tr 


The  study  will  be  based  on  a  third-order  mechanization 
error  model  of  the  vertical  loop.  In  addition,  the  analysis 
will  be  restricted  to  the  transient  response  of  the  vertical 
velocity  and  altitude  following  a  series  of  specific  man¬ 
euvers  of  the  vehicle  just  before  the  TERCOM-update .  The 
investigation  will  not  include  the  steady-state  analysis  of 
the  vertical  loop;  however,  correlation  between  the  steady- 
state  following  the  transient  behavior  will  be  analyzed. 

In  addition,  theoretical  complications  and  practical  require 
ments  will  necessitate  the  imposition  of  certain  assumptions 
1.  It  will  be  assumed  that  the  vertical  channel  can 
be  mechanized  alone.  This  means  that  the  coupling 
between  the  vertical  and  horizontal  channels  will 
be  ignored.  This  coupling  is  not  so  in  the  real 
world  environment,  however;  in  this  scope  of 
study,  it  will  not  have  a  significant  effect. 

For  a  full  scale  model,  the  coupling  between  the 
horizontal  and  vertical  channels  cannot  be  ignored. 


2.  Although  a  complete  analysis  of  the  system  requires 
that  both  transient  and  steady-state  behavior  be 
categorized,  for  the  purposes  of  analysis  in  this 
thesis,  it  will  be  assumed  that  the  vertical  loop 
of  the  INS  closely  follows  the  barometric  altitude 
in  steady-state.  This  assumption  is  true  in  the 
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practical  world,  provided  there  is  a  constant 
altitude  flight  for  a  period  greater  than  the 
time  constant  of  the  baro-inertial  vertical 
channel. 

Until  now,  the  present  day  mechanizations  of  the 
baro-aided  vertical  loop  have  been  implemented  by  concen¬ 
trating  on  the  steady-state  behavior  of  the  inertial  and 
barometric  data.  It  is  possible  that  shorter  time  con¬ 
stants  and  faster  recovery  time  may  yield  more  accurate 
instantaneous  altitude  and  velocity  at  the  expense  of 
rather  long  term  altitude  errors  such  as  those  due  to 
prolonged  descents.  This  factor  could  significantly 
improve  the  performance  of  a  vehicle  carrying  out  a 
TERCOM-type  update  following  immediately  after  a  series 
of  horizontal  and  vertical  turns.  The  present  day  mech¬ 
anizations  have  imbedded  in  them  long  time  constants  so 
that  the  INS  altitude  follows  closely  the  barometric 
altitude  in  steady-state  and  neglects  any  variations  in 
the  latter  due  to  standard  setting  or  a  scale  factor 
error.  However,  in  a  prolonged  descent/ascent,  signifi¬ 
cant  error  develops  in  the  barometric  data  with  the  INS 
closely  following  it.  Thus,  if  a  target  needs  to  be 
attacked  or  if  a  TERCOM-update  is  required  immediately 
following  a  prolonged  ascent/descent,  the  vertical  vel¬ 
ocity  and  altitude  will  be  in  significant  error. 


The  approach  will  be  to  simulate  the  vertical  channel 
and  model  the  error  propagation  in  the  presence  of  analyti¬ 
cal  models  for  the  disturbance  and  to  search  for  the  optimum 
gains  based  on  a  cost  function  which  concentrates  on  the 
TERCOM-type  measurement  update  time-frame. 

Overview 

This  thesis  is  presented  in  seven  parts.  First, 

Chapter  I  provides  a  background  and  the  necessity  for 
such  an  investigation.  Chapter  II  discusses  the  model 
selection,  including  all  uncertainties,  and  provides  a 
cost  function  in  light  of  the  mathematical  development 
present-1/.;...  Chapter  III  details  the  minimization  routine 
along  with  its  verification  and  delineates  pitfalls  and 
solutions  to  possible  problems  which  could  be  encountered 
during  this  process.  In  Chapter  IV,  the  truth  model 
and  error  state  propagation  of  the  LN-15  are  presented 
for  the  Monte  Carlo  simulation.  The  trajectory  for 
Monte  Carlo  simulation  is  also  presented  in  this  chapter. 
Chapter  V  presented  the  optimal  gains  and  also  the  vali¬ 
dation  results  from  the  Monte  Carlo  simulation.  In 
Chapter  VI,  conclusions  and  recommendations  of  this 
thesis  are  presented.  The  Appendices  contain  the  detailed 
description  of  the  instability  of  the  vertical  channel 
and  the  computer  listings  for  the  minimization  routine  and 
user  input  routines  for  the  Monte  Carlo  simulation. 


II .  Performance  Assessment  of  Vertical  Channel 

Selection  of  Model 

It  is  well  known  that,  in  the  mechanization  of  a 
pure  inertial  navigation  system,  the  calculation  of 
altitude  is  unstable  (Appendix  A  and  Ref  5) .  Several 
methods  have  been  proposed  to  stabilize  the  vertical 
channel.  Various  error  models  have  been  proposed  depend¬ 
ing  on  the  actual  application  of  the  inertial  system  in 
the  real  world  environment.  In  a  conventional  local- 
level  system,  the  stabilization  of  the  altitude  is  accom¬ 
plished  by  correcting  the  vertical  channel  integrators 
with  the  difference  between  the  inertial  system  and  alti¬ 
meter  indication  of  the  vertical  position.  Depending 
upon  the  complexity  of  the  requirement,  low-order  to 
high-order  mechanizations  are  used.  Usually,  however,  a 
third-order  mechanization  is  preferred  for  the  reasons  of 
optimum  balance  between  performance  and  mathematical 
tractability .  For  this  reason,  a  third-order  mechaniza¬ 
tion  of  the  vertical  channel  was  chosen  for  the  purposes 
of  study  in  this  thesis. 

At  present,  the  classical  third-order  mechanization 


is  in  widespread  use.  Efforts  have  been  made  toward 
improving  the  loop  gains  so  as  to  obtain  an  equitable 
balance  between  the  errors  of  the  vertical  velocity  and 


altitude.  Widnall  and  Sinha  (Ref  4)  formulated  the 
selection  of  the  three  loop  gains  in  the  baro-inertial 
vertical  channel  as  a  stochastic  optimal  control  problem, 
with  the  objective  of  minimizing  the  mean  squared  error 
of  the  indicated  vertical  velocity.  With  the  optimum 
gains  thus  obtained,  they  showed  an  improvement  of  30 
percent  over  the  classical  set  of  gains  in  a  simulated 
flight  of  an  aircraft.  A  similar  kind  of  study  was 
carried  out  in  this  thesis  with  the  objective  of  mini¬ 
mizing  the  mean  squared  error  of  the  altitude  at  TERCOM- 
update.  The  error  in  altitude  of  the  INS  at  TERCOM- 
update  is  of  far  greater  importance  than  the  error  in 
vertical  velocity.  Error  in  vertical  velocity  is  critical 
during  a  weapon  delivery  because  this  error  greatly 
affects  the  miss  distance  of  the  weapon  on  a  target. 

Since  no  weapon  delivery  is  performed  at  TERCOM-update , 
it  is  logical  to  concentrate  on  minimizing  the  altitude 
error  to  protect  against  an  incorrect  or  missed  update 
of  the  navigation  system.  In  addition  to  minimizing  the 
mean  squared  altitude  error,  a  non-stochastic  disturbance 
from  the  baro-altimeter  is  modeled  to  account  for  the 
long  term  error  introduced  during  a  descent  prior  to  TERCOM 
update.  Any  gain  selection  should  also  treat  this  non¬ 
stochastic  error  source.  By  selecting  gains  to  minimize 
altitude  error  due  to  this  disturbance  just  prior  to  the 
update,  the  vertical  channel  performance  can  be  optimized 


during  the  TERCOM-update.  The  basic  error  model  formulated 
by  Widnall  and  Sinha  was  used,  and  alterations  were  made 
according  to  the  mission  requirements,  as  will  be  shown 
in  the  succeeding  paragraphs. 

Simplified  Model 

Figure  1  shows  the  simplified  version  of  the  baro- 
inertial  vertical  channel  error  model.  The  set  of 
equations  describing  this  diagram  is: 

(1) 

(2) 

(3) 

where  6h  is  the  closed  loop  altitude  error,  6VZ  is 
the  vertical  velocity  error,  <5d  represents  the  distur¬ 
bance  input  and  is  the  variation  of  sensed  altitude  error 

A 

from  the  true  value,  and  6a  is  the  vertical  acceleration 
error  estimate  variable.  The  loop  gains  are  given  by  , 

K2  »  K3  •  an<*  9  is  magnitude  of  gravity  computed  as 
a  function  of  indicated  altiutde  and  latitude  with  R 
being  the  geocentric  radius.  The  feedback  of  6h  through 
2g/R  reflects  changes  in  gravity  with  altitude,  and  it 
can  be  recognized  as  the  cause  of  instability  in  the  unaided 


6h  =  6VZ  -  K±  <6h  -  6D) 

A 

6VZ  =  (2g/R)  6h  -  K2  (6h  -  6D)  -  6a 

A 

6a  *  K3  ( 6h  -  6D) 
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Figure  1.  Simplified  Diagram  of  Vertical  Channel 
Error  Model  with  Disturbance  Input 


inertial  vertical  channel.  Although  a  simpler  second- 
order  damping  system  would  result  from  setting  equal 

to  zero,  the  performance  advantage  of  the  more  complex 
system  is  sufficient  to  warrant  its  use. 

It  will  be  assumed  that  the  disturbance  (<$D)  acts 
on  the  vertical  channel  for  a  time  interval  At^  which 
begins  at  some  time  t^  and  terminates  at  time  t 2 
(Fig.  2) .  It  will  be  further  assumed  that  it  is  desired  to 
minimize  the  altitude  error  over  an  interval  At2  which 
begins  at  time  t^  where  t^  >  (Fig.  2) . 

In  light  of  the  above  statements,  the  cost  function 
to  be  minimized  over  At2  is 

fc4 

J(K)  =  /  ( <5h)  2  dt  ;  K  *  KlfK2,K3 

t3 

where  J (K)  is  the  performance  index  which  is  an  arbi¬ 
trary  mathematical  expression  designed  to  measure  how  well 
a  system  performs  a  particular  task.  Since  both  positive 
and  negative  values  of  the  altitude  error  are  equally 
undesirable,  the  measurement  of  the  mean  squared  error 
of  the  altitude  is  an  appropriate  means  of  indicating  how 
well  the  INS  performs  over  the  defined  time  interval. 
Another  form  of  the  cost  function  of  Eq  (4) ,  although 
mathematically  less  desirable,  would  be  to  replace  the 
square  of  the  altitude  error  by  simply  its  magnitude. 


Figure  2.  Disturbance  Profile 


Equation  (4)  is  the  basic  cost  function  and  was 
minimized  through  a  search  routine,  the  results/discussions 
of  which  are  given  in  Chapter  V. 

Addition  of  Uncertainties 

The  diagram  of  Figure  1  is  by  no  means  a  comprehensive 
depiction  of  all  the  errors  associated  with  the  vertical 
channel.  Numerous  other  error  sources  associated  with  the 
vertical  channel  must  be  modeled  to  account  for  known 
error  producing  mechanisms.  Although  the  various  error 
sources  have  been  modeled  as  white  noises  and  random 
walks,  there  are  better  models  than  just  these.  It  is 
only  for  the  sake  of  mathematical  tractability  that  the 
fir  simpler  models  are  preferred  whenever  possible.  A  compre¬ 

hensive  diagram  is  shown  in  Figure  3. 

The  feedback  path  (2g/R)  arises  from  the  gravity 
calculation  and  has  the  effect  of  destabilizing  the  altitude 
(Appendix  A) .  The  error  state  5a  is  a  random  walk  and 
it  models  the  following  (Ref  4) : 

1.  Bias  or  slowly  varying  error  in  the  vertical 
acceleration  due  to  accelerometer  bias. 

2.  Gravity  anomoly. 

3.  Error  in  Coriolis  terms. 

The  white  noise  w  ~  feeding  into  the  integrator 

a  4 

provides  the  random  walk  for  the  error  state  6a  .  The 
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2g/R 


Figure  3.  Baro-Inertial  Vertical  Channel  Error  Model 


white  noise  wa^  into  the  summing  junction  models  short 
correlation  time  accelerometer  error,  which  could  arise 
due  to  the  vertical  accelerometer  scale  factor  error  and 
input  axis  misalignment  during  the  maneuver  of  the  vehicle 
(Ref  4).  The  random  white  noise  w^  models  any  short 
correlation  time  altimeter  error  due  to  change  in  side¬ 
slip  angle  or  in  angle  of  attack  during  a  maneuver  (Ref  4) , 
The  white  noise  w^  provides  the  random  walk  for  the 
error  state  5b  which  represents  the  baro-altimeter  error 


and  is  the  sum  of  terms  as  follows: 


6b  = 


+  h  e. 


V  -  Tb  vz  +  SD 


where  e  is  the  altimeter  error  due  to  altimeter  bias, 

Po 

e,  *  is  the  altimeter  scale  factor  error,  c  represents 
hsr  sp 

the  static  pressure  measurement  error  and  is  the 

altimeter  lag  during  ascents/descents.  The  additional 
term  6D  represents  the  disturbance  input  to  the  baro- 
altimeter,  which  is  present  only  during  the  time  interval 
&t2  as  stated  previously.  It  is  assumed  here  that  the 
disturbance  input  is  uncorrelated  with  all  the  four  white 
gaussian  noises.  The  modeling  details  of  the  terms  of 
Eq  (5)  are  presented  in  Chapter  IV  under  the  Truth  Model 
heading . 
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Changed  Model/Cost  Function 


01 


In  Figure  3,  the  baro-altimeter  error  6b  is  modeled 
as  a  random  walk  with  the  white  noise  w^  driving  force 
on  this.  The  random  walk  is  the  output  of  an  integrator 
driven  by  white  gaussian  noise.  Thus 

6b(t)  =  wb2 (t)  ;  6b(tQ)  =  0  (6) 

where  the  white  noise  w^  has  zero  mean  and  covariance 
dynamics . 


*6b6b(t)  “  Q 

where  Q  is  the  strength  of  the  white  noise  driving  the 
integrator . 


w 


b2 


b 


It  can  be  seen  from  Eq  {7)  that  the  mean  squared  value 
grows  linearly  with  time  and  is  unbounded;  i.e. , 


E{ 6b2 (t) }  =  Q[t  -  tQ] 


(8) 
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Since  no  term  of  Eq  (5)  grows  unbounded  with  time,  it  is 
inappropriate  to  model  6b  as  a  random  walk.  For  this 
reason,  this  model  was  deleted  and  instead  6b  was 
modeled  as  the  output  of  a  first-order  lag  driven  by 
white  gaussian  noise  as  follows  (Ref  8)  (the  selection  of 
T  is  done  in  the  next  section) : 


This  model  produces  an  autocorrelation 


'^6b6b^'0  =  E{6b(t)6b(t  +  x)}  =  o2e 


-  I  x  I  /T 


i.e.,  of  correlation  time  T  and  mean  squared  value  a2 
(with  zero  mean)  (the  selection  of  T  and  o  is  done  in 
the  next  section) .  Thus,  the  first  order  lag  can  be 
described  as 


6b(t)  =  -  (1/T) 6b (t)  +  w ( t ) 
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r 
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'a' 


rs 


I  *» 


where  Q  is  the  strength  of  the  white  noise  and  the 
mean  squared  value  of  the  process  is 


where 


E{ 6b2 { t) }  =  =  a2 


(ID 


Q  =  2c2/T  (12) 

It  can  be  clearly  seen  from  Eq  (11)  that  the  variance  is 
constant. 

The  changed  model  of  the  baro-inertial  vertical  channel 
is  depicted  in  Figure  4,  where  the  baro-altimeter  error 
is  modeled  as  the  output  of  a  first-order  lag  driven  by 
white  gaussian  noise  w^2 

In  the  presence  of  random  inputs,  as  is  the  case 
here,  it  is  more  appropriate  to  take  the  expected  value 
of  the  cost  function  denoted  as  E{*}  herein.  Recalling 
Eq  (4),  the  cost  function  can  be  written  as 


J(K) 


E{/ 

t- 


6h2 (t)  dt) 


(13) 


Define 


6h(t)  =  6h1(t)  +  6h2(t) 


(14) 
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Figure  4.  Changed  Model  of  Baro-Inertial  Vertical  Channel 


where  Sh^t)  and  Sh2(t)  are  the  errors  in  altitude 
due  to  the  deterministic  (disturbance)  and  white  gaussian 
noises  (of  mean  zero) ,  respectively.  Taking  the  square 
of  Eq  (14)  yields 


6h2(t)  =  6h2(t)  +  26h1 (t) 6h2 (t)  +  6h2 (t) 


Rewriting  the  right  hand  side  of  Eq  (13) 


E {/  6h2 (t) dt } 

t. 


/  E{  <Sh2  (t)  }dt 


Using  Eq  (14a) , 


E {/  6h2  (t)  dt }  =  /  E{6h.2  (t)dt)dt 

i.  i.  J* 


+2  /  E(5h1(t) 5h2 (t) }dt 


+  /  E{<5h2(t))dt 

fc3 


Since  the  disturbance  is  assumed  independent  of  the  white 
gaussian  noises  (as  stated  earlier),  Eq  (16)  simplifies  to 


E {/  6h2 (t) dt } 


/  E{ 6h2 (t ) }dt 
fc3 


+  /  E{Sh2(t)}dt 

By  our  definition  of  Sh^(t)  being  the  error  due  to  a 
deterministic  input,  the  expectation  on  6h1(t)  can  be 
removed 


J  (K) 


/  6h2 (t) dt  + 

*3 


/  E{ 6h2 (t) }dt 
fc3 


(17) 


(18) 


If  the  mean  square  value  of  the  altitude  error  due  to 
zero-mean  stochastic  inputs  remains  constant  over  the 
TERCOM-update  interval,  i.e.,  At2  (see  Fig.  2),  then 
Eq  (18)  can  be  finally  written  as 

fc4 

J  (K)  =*  /  6h2  (t)  dt  +  P6h  At  (19) 

t3 

where  is  the  covariance  (or  mean  squared  value  since 

the  stochastic  inputs  are  zero  mean)  of  the  altitude  error 
(stochastic  inputs)  over  the  interval  At  =  (t^  -  t^) 

The  first  term  on  the  right  hand  side  of  Eq  (19)  is  due  to 
the  deterministic  input  (as  in  Eq  (4) )  and  the  second  one 
—  is  due  to  the  zero  mean  stochastic  inputs.  Equation  (19) 
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can  be  written  as  the  sum  of  two  cost  functions 


fl 


J(K)  =  JL(h)  +  J2(h) 


(19a) 


where 


Jl(K) 


/  6h? (t) dt 

fc3 


and 

j2.(K)  =  P5h  At 


(19b) 


(19c) 


Thus,  the  minimization  of  the  cost  function  of  Eq  (19) 
will  lead  to  the  minimization  of  the  altitude  error  and 
the  disturbance  over  a  time  interval  At2  .  The  results 
of  minimization  of  Eq  (19)  and  discussions  are  presented 
in  Chapter  V  under  the  New  Cost  Function  heading. 

Mathematical  Development 

As  stated  in  the  previous  section,  it  is  required  to 
compute  the  mean  squared  altitude  error  as  a  function  of 
input  noise  spectral  densities  and  loop  gains.  This 
expression  is  required  in  Eq  (19)  to  compute  the  overall 
cost  function. 

It  is  useful  to  express  the  power  spectral  density 
of  a  wide-sense  stationary  output  of  a  system  directly 
in  terms  of  the  power  spectral  density  of  the  input  and 
the  description  of  the  system  itself  (Ref  8) . 
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G  ( s ) 


(s) 


^zz  ^ 


From  above,  for  a  system  of  transfer  function  G(s)  and 
input  and  output  power  spectral  densities  of  tj7nn(s)  and 
tj7zz(s),  respectively,  it  is  true  that  (Ref  8) 


^2Z  (w)  =  G  (-oj)  G  (cj)  ’^nn  (u) 


(20 


If  input  is  a  whLte  gaussian  noise  of  strength  Q  for 
all  oj  ,  then  Eq  (20)  becomes 


^Zz(oj)  =  G(-w)G(w)Q 


(21 


Similar  to  the  lines  of  the  development  of  Eq  (21) ,  the 
power  spectral  density  of  the  altitude  error  can  be 
computed  from  the  following  equation: 

4 

^6h(u)  =  E  G±  ( j  u))  Gi(-ju)  Qi  (22 


or 


WS) 


4 

l  G. (S)  G. (-S)  Q. 

•  i  1  X  X 


(22 


where  G.  is  the  transfer  function  from  each  of  the 
1 

independent  white  noise  sources  to  the  output  (Fig.  3) 
and  Qi  is  the  strength  of  the  individual  white  noises 
of  all  the  four  sources.  The  mean  squared  value  of  the 
altitude  error  due  to  stochastic  inputs  is  then  the  inte 
gral  of  the  power  spectral  density 

<5h2»2  -  nr[^Gi(s»  Gi(-S>  ds 


To  calculate  the  mean  squared  altitude  error,  we  first 
need  to  calculate  the  trnasfer  function  from  each  of  the 
individual  white  noise  sources  to  the  output.  This 
simple  expression  results  from  the  assumption  that  the 
white  noise  sources  are  uncorrelated  and  independent, 
since  each  error  is  from  a  different  source  and  there 
is  not  reason  to  believe  that  they  are  correlated. 

Figure  5  shows  the  various  transfer  function  blocks 

between  the  white  noise  w  .  and  5 h 

al 

The  overall  transfer  function,  for  the  first  noise 


Gi  ( s) 


_ S _ 

K,  +  (K,  -  C)S  +  S  (S  +  K.) 


Figure  5.  Transfer  Function  Block  from  White  Noise  w  to  Output 
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G1(S) 


6  h 

Wal  S  +  i^S  +  (K2  -  C)S  +  K3 


With  Q  ^  being  the  strength  of  white  gaussian  noise 
w  .  and  using  Eq  (23)  for  i=l  ,  we  get 

a.  J. 

00 

'  55J  [.  G1(S)G1(-S)  dS 

Using  the  table  of  integrals  (Ref  9) ,  we  get 


(6h  )2  .  { - i - }  Q 

1  2  K,  (K_  -  C)  -  K, 


al 


1  " 

»  . 

* 


m 


For  the  second  noise  source  w^2  >  an  additional  integrator 

to  Figure  5  is  needed,  which  is  shown  in  Figure  6.  From 
Figure  6,  the  transfer  function  from  white  gaussian  noise 
wa2  strength  Qa2  to  the  output  is 

G  (S)  =  =  - - - 

Wa2  S[S  +  KlS  +  (K2  -  C)S  +  K3] 

As  before,  using  Eq  (23)  for  i=2  with  Q  0  being  the 
strength  of  w&2  , 

Q  *  °° 

<«h2>*  =  2lf  /’-  G2<S,G2(-S>  dS 
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(27) 


(28) 


(29) 


Final  Transfer  Function  for  Noise  Source  w  „  to  Output 


Using  the  table  of  integrals  (Ref  9) 


in 

i 

i 

► 

9  ,T 

P 


r  ; 

i  '' 
i 

> ‘  • 


(6h.)  2 
z  2 


{ - 

2[K1K3 (K2 


K, 


-  C) 


} 


To  calculate  the  transfer  function  for  noise  source  w 


bl 


to  output,  v/e  proceed  in  a  similar  fashion  as  before. 
Figure  7  shows  the  block  diagram  from  noise  source 
to  output.  The  overall  transfer  function  for  the  third 


noise  source  w 


bl 


is 


G,  (S)  = 

3  w. 


6h 


K,  S2  +  K-S  +  K. 


bl 


S3  +  K1S2  +  (K2  -  C)S  +  K3 


with  being  the  strength  of  wfel  and  using  Eq  (23) 

for  i=3 


(6h9) 2 
*  3 


/  G3(S)G3(-S)  ds 

~j°o  J 


Using  the  integral  tables  (Ref  9) 


(6h_) 2 
1  3 


K2 (K-  -  C)  +  K2  -  K.K- 

- 2 - - — — }  Q.. 

2(K1(K2  -  C)  -  K3] 


Lastly,  for  the  noise  source  Wja2  ,  an  addition  of  the 
first-order  lag  to  Figure  7  gives  the  desired  transfer 
function  and  is  depicted  in  Figure  8.  From  Figure  8, 


(30) 


(31) 


(32) 


(33) 


Figure  7.  Transfer  Function  Block  from  Noise  Source  w,  ,  to  Output 


Figure  8.  Final  Transfer  Function  Block  from  Noise  Source  w,  „  to  Output 


(34) 


6h  _  KXS2  +  K2S  +  K3 

Wb2  (S  +  a) [S3  +  KlS2  +  (K2  -  C)S  +  K3] 

K1S2  +  K2S  +  K3 
S4+(a+K1)S3  +  (aK1+K2-C)S2  +  (aK2-aC+K3)S  +  aK3 

(35) 


With  0^2  being  the  strength  of  w^2  and  using  Eq  (23)  for 
i=4 


( 6h_ )  2 
^  4 


Q 


b2 


2nj 


.  CO 


G4(S)G4(-S)  dS 


(36) 


Using  the  integral  tables  (Ref  9) ,  the  solution  to 
Eq  (36)  is 


( <5h9 )  2 
^  4 


ra2[K2 (K2-C) -K1K3+K2]  +  aKxK2 

+k3[k1(k2-c)  -  k3] 


lh2 


2a 


l"a3(K1  (K2-C)  -K3}+  a2(K1  (KL  (K2~C)  -K3)  } 

+  a(K1K2(K2-C)  -  K3 (K2-C) -CK1 (K2-C) } 

L-  +  k3Ck1(k2-c)  -  k3> 


=T  (37) 


Combining  Eqs  (27),  (30),  (33)  and  (37)  and  substituting 
in  Eq  (23) ,  the  mean  squared  altitude  error  becomes 
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fl 


(Sh2)2 


Q 


al 


(Kl)Qa2 


2[K1(K2-C)  -  K3]  2[K1K3(K2-C)  -  K2] 


(K2(K2-C)  +  K2  -  K1I<3)  Qbl 
2[K1(K2-C)  -  K3] 


+ 


a2[K2(K2-C)  -  KxK3  +  K2]  +  al^K 


k3[k1(k2-c)  -  k3] 


2a 


a3{K1 (K2-C) -K3}+a2 (K1 (K^K^-C) -K3) } 
+  a{K1K2 (K2-C) -K3 (K2-C) -CKl (K2-C) } 

+  k3{k1(k2-c)-k3> 


Equation  (38)  is  then  the  required  mean-squared  altitude 
error  for  use  in  Eq  (19)  to  make  the  complete  cost  function. 

One  more  aspect  still  remains  untouched.  Appropriate 
values  for  the  strength  of  the  white  gaussian  noises  and 
the  correlation  parameter  of  the  first-order  lag  are 
required.  It  is  rather  difficult  to  suggest  values  which 
provide  a  true  depiction  of  the  real  world  environment. 
However,  without  these  values,  further  progress  will  not 
be  possible.  Table  1  shows  the  nominal  values  of  the 


(38) 
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noise  spectral  densities  and  correlation  parameter  that 
have  been  selected.  Three  out  of  the  five  values  were 
selected  based  on  the  reasoning  of  Widnall  and  Sinha 
(Ref  4) .  These,  and  the  reasoning  for  the  remaining  two, 
are  described  in  the  following  section. 


As  suggested  by  Widnall  and  Sinha  (Ref  4) ,  a  typical 
root  mean  square  (rms)  amplitude,  for  a  short  correlation 
time  acceleration  error,  is  about  200  ug.  This  figure 
is  appropriate  assuming  a  horizontal  maneuver  of  duration 
of  60  seconds.  This  error  could  be  caused  by 

(1)  A  200  yrad  misalignment  of  the  input  axes 
caused  by  the  vertical  accelerometer,  and 

(2)  a  horizontal  maneuver  acceleration  of  about 
one  g. 

Assuming  a  repeated  random  maneuver,  the  area  of  the 
acceleration  error  autocorrelation  is  (as  derived  in 
Ref  4) 

Qal  =  2.4  x  10-4  m2/sec3 

The  area  of  the  autocorrelation  is  the  low  frequency 
value  of  the  power  spectral  density.  For  a  white  gaussian 
noise  whose  autocorrelation  is  the  dirac  delta  function 


with  area  Qa^  ,  the  spectral  density  applies  at  all 
frequencies  (Ref  4) .  Since  we  are  interested  in  the 
lower  frequencies  of  the  short-correlation  acceleration 
error,  therefore,  the  low  frequency  density  of  Eq  (39)  is 
used  for  the  spectral  density  of  the  white  noise  for  all 
frequencies. 

The  acceleration  error  6a  (Fig.  4)  models  the 
inertial  vertical  acceleration  error  and  it  is  caused  by 
reasons  already  outlined  earlier  in  this  chapter.  For  an 
assumed  period  of  1000  seconds,  if  the  rms  value  of  the 
accelerometer  bias  is  expected  to  shift  100  pg  approxi¬ 
mately,  then  the  strength  of  the  white  noise  wa2  as 
derived  in  Ref  4  is 


=  1.0  x  10  9  m2/sec5 


For  a  short  correlation  time  altimeter  error,  it  is 
assumed  that  an  rms  error  of  10m  may  be  present  in  the 
baro-altitude  with  a  correlation  time  of  one  second 
(Ref  4) .  Thus  the  strength  of  the  white  gaussian  noise 


as  derived  in  Ref  4  is 


=  100  m2sec 


As  explained  earlier  in  this  chapter,  the  white  noise 
Wj32  models  the  error  state  6b  (Fig.  4)  which  represents 
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the  baro-altimeter  error  which  is  the  sum  of  many  terms 
(Eq  (5)).  For  the  kind  of  trajectory  (see  Chapter  IV) 
and  for  the  minimization  of  the  altitude  error  at  the 
TERCOM-update,  the  altimeter  bias  or  the  standard  setting 
error  is  of  primary  concern.  It  is  assumed  that  the 
vehicle  has  been  in  flight  for  a  sufficiently  long  time 
over  a  great  distance  before  the  TERCOM-update  and  the 
effect  of  the  standard  setting  error  is  predominant.  If 
the  altimeter  bias  is  represented  by  e  (consistent 
with  Eq  (5)),  then  this  error  can  be  modeled  as  a  first- 
order  Markov  process  given  by  (Ref  4) 


e  =  -a  e  „  +  w.  ~ 
po  po  b2 


a  =  V/d. 


(42a) 


Qb2  2  a  aalt 2 


(42b) 


where  d  is  the  correlation  distance  of  the  weather 

system,  a  .  is  the  standard  deviation  of  the  variation 

cl  1C 

in  altitude  of  a  constant  pressure  surface,  Q ^  is  the 
power  spectral  density  of  the  white  gaussian  noise  w^  • 
and  V  is  the  vehicle  speed.  For  a  vehicle  speed  of 
600  miles/hr  (more  appropriate  for  a  missile) ,  correlation 
distance  (d^^.)  of  250  nautical  miles  and  a  one-sigma 
value  (aa]_t)  of  500  feet  (Ref  4)  ,  the  strength  of  the 


nr 


white  noise  and  the  value  of  the  correlation  parameter 
become 


a 


600  x  52S 0  -1 

-  sec 

3600  x  250  x  6080  20 


or 


a  =  5.793  x  10_4/sec 


and 


Qb2  =  (2)  (5.79  x  10-4)  (152.4  m) 2 


or 


Qb2  =  26.91  m2sec 

Equations  (43)  and  (43a)  are  valid  only  for  the 
constant  velocity  of  600  mi/hr,  and  the  value  of  the 
correlation  parameter  (a)  will  change  with  change  in 
velocity  of  the  vehicle,  as  for  example  during  descents 
or  ascents. 

As  stated  previously.  Table  1  shows  the  values  of  the 
spectral  densities  and  correlation  parameter  for  use  in 
Eq  (38)  . 


(43) 


(43a) 


(44) 


(44a) 
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TABLE  II-l 

Nominal  Values  of  Noise  Spectral  Densities 

and  Correlation  Parameter 

White  Noise  for 

NOise 

Density/ 

Correlation 

Parameter 

Value 

Short  Correlation  Time 

Acceleration  Error 

Qal 

2.4xl0-4  m2sec~3 

Acceleration  Error 

Random  Walk 

Qa2 

-  9  o  -  5 

1.0x10  mzsec 

Short  Correlation  Time 

Altimeter  Error 

Qbl 

100  m2sec 

Altimeter  Error 

First-Order  Lag 

Qb2 

26.91  m2sec-1 

Correlation  Time 

for  First-Order  Lag 

a 

5.793xl0-4  sec-1 
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III.  Program  for  Minimization  of  Cost 
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Selection  and  Development  of  Routine 

The  previous  sections  dealt  with  the  development  of 
the  cost  function  as  given  in  Eq  (19) .  The  parameters 
which  need  to  be  optimized  are  the  three  loop  gains  (K-^  , 
K2  ,  K.j)  of  the  vertical  channel.  For  convenience's 
sake,  the  cost  function  of  Eq  (19)  is  reproduced  as 

j(k1,k2,k3)  =  j1(k1,k2,k3)  +  j2(k1,k2,i<3) 

J(K1,K2,K3)  =■  /  (6h1)2dt  +  (t4  -  t3)  (P6h) 

t3 


where  the  second-half  portion  (i.e.,  the  covariance  P^ 
or  (<5h2)  2  since  the  mean  of  stochastic  inputs  is  zero) 
of  the  right  hand  side  is  given  by  Eq  (38)  and  the  first 
one  by  Eqs  (1) ,  (2) ,  and  (3) .  To  achieve  the  total  cost, 
Eqs  (1)  through  (3)  need  to  be  integrated  for  each  interval 
of  time  to  which  Eq  (38)  is  added.  Thus,  an  integration 
package  is  required  along  with  the  search  routine. 

The  integration  and  search  routines  selected  for 
this  thesis  were  DGEAR  and  ZXMIN,  respectively,  both  of 
which  reside  in  the  IMSL  Library  (reference) .  The  DGEAR 
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(45) 


(45a) 


routine  finds  approximations  to  the  solution  of  a  system 
of  first  order  ordinary  differential  equations  with  initial 
conditions.  The  basic  method  used  for  the  solution  is 
of  implicit  linear  multistep  type.  This  routine  is  very 
useful  in  solving  the  stiff  differential  equations  which 
were  encountered  during  the  course  of  this  thesis  (small 
step  sizes  were  taken  by  the  integration  routine  to  achieve 
reasonable  accuracy  for  extremely  large  values  of  the 
gain  K1  ) .  References  10  and  11  can  be  consulted  for 
more  details.  The  search  routine  ZXMIN  is  based  on  the 
Harwell  library  routine  VA10A  and  utilizes  the  quasi-Newton 
method  to  find  the  minimum  of  a  function.  The  search 
routine  ZXMIN  was  selected  because  it  requires  no  explicit 
gradient  information  from  the  user  (it  internally  computes 
the  gradient  if  not  available) .  Reference  11  and  the  IMSL 
package  can  be  consulted  for  additional  information. 

A  simple  flow  chart  of  the  computer  program  is  shown 
in  Figure  9.  Estimates  of  the  three  loop  gains  are  fed 
into  the  search  routine  which  outputs  the  values  of  the 
cost  and  the  three  loop  gains.  The  search  routine  iter¬ 
atively  estimates  the  values  of  the  loop  gains  until  a 
minimized  cost  is  obtained.  The  convergence  condition  is 
satisfied  if,  on  two  successive  iterations,  the  parameter 
estimates  (i.e.,  K^,  K^)  agree  component  by  component 

to  the  number  of  significant  digits  specified  (3  to  5 


SET  INITIAL  VALUES 


Figure  9.  Flow  Diagram  of  Search  Routine 
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significant  digits  for  this  study) .  A  sample  listing 
of  the  program  can  be  found  in  Appendix  B. 


Validity  Check  of  Program 

The  first  step  after  the  development  of  the  minimiza¬ 
tion  program  was  to  validate  it  through  an  earlier  published 
result.  The  results  obtained  by  Widnall  and  Sinha  (Ref  4) 
were  selected  for  this  comparison. 

To  do  this,  the  cost  function  becomes 


J(v)  =  ( 6v)  2  (46) 


where  (5v) 2  is  given  by  an  equation  developed  along  the 
lines  similar  to  Eq  (38)  for  the  case  of  the  mean  squared 

value  of  the  vertical  velocity  error  as  done  by  Widnall 

$ 

and  Sinha  (Ref  4) . 

Using  the  values  for  the  strength  of  the  white  noise 
as  given  by  Widnall  and  Sinha  (Ref  4) ,  the  results  obtained 
using  the  method  developed  above  for  the  three  loop  gains 
were  exactly  the  same  as  those  obtained  by  them,  thus 
validating  the  minimization  program. 

Scaling  and  Techniques  Used 

It  is  well  known  that  one  of  the  greatest  pitfalls 
of  a  computer  search  routine  is  that  the  routine  is  liable 
to  converge  toward  a  local  minimum,  whereas  what  is  needed 
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is  the  global  minimum  of  the  function.  Unless  the 
function  to  be  minimized  is  well  defined  (in  which  case 
the  local  minimum,  different  from  a  global  minimum,  does 
not  exist) ,  the  results  obtained  from  a  computer  search 
are  often  questionable.  To  overcome  this  problem,  it  is 
advisable  to  have  many  sets  of  starting  points  for  the 
input  variables. 

For  this  thesis,  the  input  variables,  as  stated 
earlier,  are  the  three  loop  gains  (K^,  I<2 ,  )  of  the 

vertical  channel.  Dr.  Widnall  and  Sinha  (Ref  4)  dis¬ 
covered  through  dimensional  analysis  a  correct  expression 
that  also  gives  an  approximate  value  for  the  gain  : 


where  and  are  the  strength  of  the  white  noises 

associated  with  wfal  and  w^  ,  respectively  (see  Fig.  3) 
With  the  approximate  value  of  known,  it  was  relatively 

easy  to  implement  the  search  routine.  Nonetheless,  differ¬ 
ent  starting  points  for  gains  K2  and  were  tried  to 

ensure  a  global  minimum.  Since  this  routine  also  utilized 
an  integration  package,  it  was  well  worth  the  effort  to 

_»  O 

keep  a  tight  control  on  the  tolerance  (  <  10  )  in  the 

integration. 


Another  problem  encountered  was  the  fact  that  the 
values  of  the  three  loop  gains  differ  from  each  other 
by  a  large  magnitude.  Consequently,  it  was  necessary  to 
scale  the  three  gains  to  the  same  level  before  they  were 
fed  into  the  search  routine.  This  procedure  ensured  that 
the  same  number  of  significant  digits  was  obtained  in  the 
final  values  of  t^e  loop  gains  and  it  also  simplified  the 
job  of  the  optimization  algorithm. 

With  the  cost  developed  in  Chapters  I  and  II,  and  the 
search  routine  developed  above,  values  for  the  optimum 
loop  gains  were  found  (see  Chapter  V) .  It  was  now  neces¬ 
sary  to  validate  these  optimum  gains  through  a  simulated 
flight  of  a  vehicle  which  is  presented  in  the  following 
chapter. 


IV.  Error  State  Propagation  and  Simulation 


rt 


The  Truth  Model 

A  "truth  model"  is  the  analytic  designer's  best 
description  of  the  real  world  behavior  of  the  INS.  In 
this  section,  a  50  state  system  error  model  (or  truth 
model) ,  which  is  needed  for  Monte  Carlo  study  of  optimal 
gains  (covariance  analysis  program  was  not  used  due  to 
non-zero  mean  disturbance  since  one  of  the  requirements 
of  covariance  analysis  is  that  all  inputs  should  be  of 
zero  mean) ,  is  presented  in  the  form  of  a  stochastic 
linear  vector  differential  equation  as  shown  in  Eq  (48) . 

x(t)  =  F  (t )  x  (t)  +  G  (t)  w(t) 

where 

x(t)  is  the  50  dimensional  state  vector, 

F (t )  is  the  (50x50)  error  propagation  matrix, 

w(t)  is  a  (10x1)  vector  of  white  noise  forcing 
functions,  and 

G (t )  is  a  (50x10)  input  matrix. 

The  error  model  of  50  state  variables  documented  in  this 
thesis  is  the  Litton  LN-15  navigation  system  (Ref  7)  with 


(48) 
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the  platform  oriented  in  East-North-Up  (ENU)  local-level 
frame.  The  50  state  variables  are  presented  in  Table  IV-1. 
Variables  1  to  9  are  the  basic  position  velocity  and 
attitude  variables.  Variable  10  is  the  additional  inte¬ 
gration  in  the  altitude  channel  mechanization.  Variables 
11  to  43  are  the  gyro  and  accelerometer  sources  of  error. 
Variables  44  to  50  are  the  altimeter  errors  and  gravity 
disturbances.  Variables  11  to  50  are  modeled  as  random 
constants,  random  walks  and  first-order  Markov  processes. 
The  models  are  briefly  summarized  in  this  thesis.  Details 
are  available  in  Reference  8 . 

A  random  constant  is  modeled  as  the  output  of  an 
integrator  with  zero  input  and  an  initial  condition  which 
has  a  zero  mean  (could  be  non-zero  mean)  and  a  variance 

P  .  The  model  is  suitable  for  an  instrument  bias  that 
o 

changes  each  time  the  instrument  is  turned  on,  but  remains 
constant  while  the  instrument  is  on. 

The  random  walk  model  is  the  output  of  an  integrator 
driven  by  a  zero  mean  white  gaussian  noise.  The  defining 
equations  are 

x  (t)  =  w  ( t )  x(tQ)  =  0 

E  {  w  ( t )  }  =  0 

E{w  (t)  w  (t+T )  }  =  Q(t)5(i) 


(49) 

(50) 

(51) 
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TAELE  IV-1 


Error  Model  State  Variables 


Basic  Inertial  Navigation  Errors 


1. 

6A 

Error 

in 

east  longitude 

2. 

6L 

Error 

in 

north  latitude 

3. 

6h 

Error 

in 

altitude 

4. 

6Ve 

Error 

in 

east  velocity 

5. 

6V 

n 

Error 

in 

north  velocity 

6. 

5V 

z 

Error 

in 

vertical  velocity 

7. 

e 

e 

East  attitude  error 

8. 

e 

n 

North 

attitude  error 

9. 

ez 

Vertical 

attitude  error 

Vertical  Channel  Error  Variable 

/■s 

10.  <$a  Vertical  acceleration  error  variable  in 

altitude  channel 

G-Insensitive  Gyro  Drifts 


11. 

DXf 

x-gyro 

drift 

rate 

12. 

°Yf 

y-gyro 

drift 

rate 

13. 

DZf 

z-gyro 

drift 

rate 

TABLE  IV- 1  (Cont'd) 


G-Sensitive  Gyro  Drift 


14. 

DX 

X 

x-gyro 

input 

axis  g-sensitivity 

15. 

DX 

y 

y-gyro 

spin 

axis  g-sensitivity 

16. 

Dy 

1  X 

y-gyro 

spin 

axis  g-sensitivity 

17. 

Dy 

y 

y-gyro 

input 

axis  g-sensitivity 

18. 

DZ 

y 

z-gyro 

spin 

axis  g-sensitivity 

19. 

DZ 

z 

z-gyro 

input 

axis  g-sensitivity 

G2-Sensitive  Gyro  Drift  Coefficients 


20. 

DX 

xy 

x-gyro 

spin 

input 

g 2 -sensitivity 

21. 

DY 

xy 

y-gyro 

spin 

input 

g2-sensitivity 

• 

C*4 

(N 

DZ 

yz 

z-gyro 

spin 

input 

g2-sensitivity 

Gyro  Scale  Factor  Errors 


23. 

GSF 

X 

x-gyro 

scale 

factor 

error 

24. 

GSF 

y 

y-gyro 

scale 

factor 

error 

25. 

GSFZ 

z-gyro 

scale 

factor 

error 

Gyro  Input  Axis  Misalignments 


26. 

XG 

y 

x-gyro 

input 

axis 

misalignment 

about 

y 

27. 

XGz 

x-gyro 

input 

axis 

misalignment 

about 

z 

28. 

YGx 

y-gyro 

input 

axis 

misalignment 

about 

X 

29. 

YGz 

y-gyro 

input 

axis 

misalignment 

about 

z 

30. 

ZG 

X 

z-gyro 

input 

axis 

misalignment 

about 

X 

31. 

ZG 

\T 

z-gyro 

input 

axis 

misalignment 

about 

y 
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TABLE  IV-1  (Cont'd) 


Accelerometer  Biases 


32. 

ABx 

x-accelerometer 

bias 

33. 

AB 

y 

y-accelerometer 

bias 

34. 

AB 

z 

z- accelerometer 

bias 

Accelerometer 

Scale  Factor  Errors 

35. 

ASF 

X 

x-accelerometer 

scale 

factor  error 

36. 

ASF 

y 

y-accelerometer 

scale 

factor  error 

37. 

ASF 

z 

z -accelerometer 

scale 

factor  error 

Accelerometer 

Input  Axis  Misalignment 

38. 

XA 

x-accelerometer 

input 

axis 

misalignment 

y 

about  y 

39. 

XAZ 

x-accelerometer 

input 

axis 

misalignment 

about  z 

40. 

YA 

y-accelerometer 

input 

axis 

misalignment 

X 

about  x 

41. 

YAz 

y-accelerometer  input  axis  misalignment 

about  z 

42. 

ZA 

z -accelerometer 

input 

axis 

misalignment 

X 

about  x 

43. 

ZA 

V 

z -accelerometer 

input 

axis 

misalignment 

about  y 


Barometric  Altimeter  Error 


44 


po 


Error  due  to  variation  in  altitude  of 
a  constant  pressure  surface 
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TABLE  IV-1  (Cont'd) 


Gravity  Uncertainties 


45.  Sg, 


East  deflection  of  gravity 


46.  6gn  North  deflection  of  gravity 


47.  6g, 


Gravity  anomaly 


*Additional  Baro-Inertial  Altimeter  Errors 


48 .  e. 


Scale  factor  error 


49.  C 


Coefficient  of  static  pressure 
measurements 


50.  Tl 


Altimeter  lag 


*These  states  were  grouped  separately  from  state  44 
to  conform  to  those  given  in  Reference  7. 
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A  first  order  Markov  model  is  the  output  of  a  first 
order  lag  driven  by  a  zero  mean  white  gaussian  noise  of 
strength  Q  .  The  equations  are 

x  (t)  =  -  x  (t)  +  w(t)  ( 

E [x2 (t ) ]  =  QT/2  ( 

where  T  is  correlation  time.  The  first  order  Markov 
model  is  a  useful  shaping  filter,  providing  adequate  approx¬ 
imation  to  a  wide  variety  of  empirically  observed  band- 
limited  (wide  or  narrow  band)  noises. 

The  error  propagation  matrix  of  the  vector  differential 
equation  governing  the  50  state  variables  is  presented  in 
partitions  in  Figures  10,  11,  12,  and  13.  Figure  10  presents 
the  upper  (9x9)  fundamental  matrix  (Pinson  error  model)  of 
INS  general  error  differential  equations.  Figure  11  shows 
those  elements  that  must  be  added  to  the  elements  of  the 
general  9x9  error  propagation  matrix  of  Figure  10  to  obtain 
the  partition  of  position,  velocity,  attitude  and  vertical- 
acceleration  error  state  variables.  Figure  12  presents  the 
non-zero  elements  of  the  gyro-error  columns  of  the  error 
propagation  matrix.  Figure  13  presents  the  non-zero  elements 
of  the  accelerometer,  gravity  disturbance  and  altimeter 
columns  of  the  fundamental  matrix.  Notation  used  in  the 
above  mentioned  figures  is  defined  in  Table  IV-2.  Error 
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igure  12.  Partition  of  Gyro-Error  State  Variables 
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Partitions  of  F-Matrix  of  Accelerometer 
Altimeter  and  Gravity  Disturbance  State 


TABLE  IV-2 


Notation  Used  in  Figure  10  to  Figure  13 


Latitude  of  Vehicle 


Earth  rotation  rate 


Radius  of  Earth 


Gravity  vector  magnitude 


V  ,  V  ,  V 
e  n  z 

f  ,  f  ,  f 
e'  n  z 


=  ficosL 
=  ftsinL 


Vehicle  velocity  with  respect  to  earth 
Components  of  specific  force 


Components  of  earth  rate 


VR 

Ve/R 


tanL 


Components  of  angular  velocity  of  ENU 
frame  with  respect  to  earth 


*  pe 

un  -  Pn+^n  Components  of  angular  velocity  of  ENU 

1  ^  frame  with  respect  to  inertial  space 

-  p 

z  z  z 


Kz  -  VR 


2 <a  v  +a  v  )  +  p  v  /cos2l 
n  n  z  z  ^n  n 


p  o  +  p  K 
Kz  e  n  z 


■p  tanL  -  K 
Ke  z 


-2ft  V  -  p  V  /cos ZL 


TABLE  IV- 2  (Cont'd) 


pnpz  "  PeKz 
2g/R  -  (p*  _  p») 

to  +  p  tanL 
n  z 

Vertical  channel  loop  gains 
Height  of  the  vehicle  over  earth 
cosa,  cosine  of  wander  angle 

sina,  sine  of  wander  angle 

sina,  sine  of  wander  angle 

cosa,  cosine  of  wander  angle 

Components  of  specific  force 

f  cos  a 
x 

f  sin  a 
f  cos  a 
f  sin  a 

y 

-f  sin  a 

f  cos  a 

-f  sin  a 

y 

f  cos  a 

y 

f  f  cos  a 
x  y 


f  f  sin  a 


TABLE  IV-2  (Cont'd) 


r? 


F721 

= 

-f  f  sin  a 
x  y 

F821 

= 

f  f  cos  a 
x  y 

F922 

= 

f  f 

y  z 

^x'  uy 

” “ 

Components  of  angular  velocity  of  ENU 
frame  with  respect  to  inertial  space 
along  LN-15  x,y  axes 

^z 

= 

Up  component  of  earth  rate 

dalt 

= 

Correlation  distance  of  altimeter  error 

d  ,d  ,d 
ge  gn  gz 

= 

Correlation  distances  of  gravity 
deflections  and  anomaly 

alt 

= 

la  amplitude  of  altimeter  error  e^ 

a  ,  a  ,  a 
ge  gn  gz 

= 

la  amplitude  of  gravity  disturbances 

F723 

= 

u  c°s  a 

F823 

= 

s  m  ct 

F724 

= 

-coy  sin  a 

F824 

= 

a)  cos  a 

F726 

= 

ftz  cos  a 

F826 

- 

ftz  sin  a 

F727 

= 

-w  cos  a 

F827 

= 

-w  sin  a 

F72  8 

= 

ft  sin  a 

z 

F828 

= 

-ft  cos  a 
z 
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p 

729 

= 

-w  sin  a 

X 

p 

£  829 

= 

wx  cos  a 

V 

= 

ground  speed  of  vehicle 

F 

349 

= 

KjV2 

F 

£  350 

= 

-K.V 

1  2 

F432 

= 

cos  Ot 

F532 

= 

sin  a 

F433 

= 

-sin  a 

F533 

= 

cos  a 

F435 

= 

f  cos  a 

A 

F535 

= 

f  sin  a 

F436 

- 

-f  sin  a 

y 

F536 

- 

f  cos  a 

F438 

= 

-f  cos  a 
z 

F538 

s 

-f  sin  a 
z 

F439 

= 

f  cos  a 

y 

P 

539 

= 

f  sin  a 

F440 

= 

-f  sin  a 
z 

F540 

= 

f  cos  a 
z 

F441 

f  sin  a 

X 

59 


TABLE  IV- 3 


Error  Source  Initial  Values  and  Statistics 


Random  Walks  x  =  W 

State 

Variable 

Error  Source 

Initial 

Value 

Noise 

Spectral 

Density 

11,  12 

x  and  y  gyro  drifts 

0 . 5xl0~3  °/hr 

(0°/hr) 2 /hr 

13 

z  gyro  drift 

0.7x10“ 3  °/hr 

(0°/hr) 2/hr 

32,  33,  34 

x,  y  and  z  accelerometer  biases 

25  yg 

(0  yg)2/hr 

First  Order  Markov  Processes 


x  =  -8x  +  W  ;  N 

w  =  23°2 

State 

Variable 

Error  Source 

Initial 

Value 

Noise 

Spectral 

Density 

Correlatior 

Parameter 

44 

Baro  Altimeter  Bias 

50  ft 

500  ft2/nm 

250  nm 

45 

East  Gravity  Deflection 

26  Mg 

140  ug2/nm 

10  nm 

46 

North  Gravity  Deflection 

17  Mg 

5  8  iiy  2/nm 

10  nm 

47 

Gravity  Anomaly 

35  ug 

41  ug2/nm 

60  nm 

(Continued  on  next  page) 
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TABLE  IV- 3  (Cont’d) 


• 

Random  Constants  x  =  0 

State 

Variable 

Error  Source 

Initial 

Value 

14 

to 

19 

G-sensitive  gyro  drift 
coefficients 

0.5xl0-3  °/hr/g 

20 

to 

22 

G2 -sensitive  gyro  drift 
coefficients 

0  °/hr/g2 

23 

to 

25 

x,  y  and  z  gyro  scale  factor 

5  ppm 

26 

to 

31 

Gyro  input  axis  misalignment 

2.5  arc  sec 

35 

to 

37 

Accelerometer  scale  factor 
errors 

25  ppm 

38 

to 

43 

Accelerometer  input  axis 
misalignments 

5.51  arc  sec 

48 

Altimeter  scale  factor 

0.003 

49 

Static  Pressure  Measurement 
Error 

0 . 1540xl0~ 3  ft 
(ft/sec) 2 

50 

Altimeter  lag 

0.25  sec 

source  initial  values  and  statistics  are  summarized  in 
Table  IV-3 . 

Not  included  in  the  above  is  the  effect  of  the  dis¬ 
turbance  to  the  baro-inertial  vertical  channel.  For  this, 
it  is  necessary  to  show  how  each  error  associated  with  the 
vertical  channel  is  modeled. 

Recalling  from  Chapter  I,  the  closed  loop  altitude 
error  6h  ,  the  vertical  velocity  error  6vz  and  vertical 

A 

acceleration  6a  estimate  are 


5h  =  6v  -  K  ( 6 h  -  6h.  ) 

Z  X  D 

5vz  =  (^)6h  -  K2(6h  -  6hfa)  -  6a 

x 

6a  =  K3(6h  -  6hb) 


where  6hfa  is  the  baro-altimeter  error  and  is  given  by 


e  _  +  he.  e 
oo  hsf 


+ 


T,  v  +  6D 
b  z 


where 


vehicle  altitude 
vehicle  speed 
vertical  velocity 


(55) 

(56) 

(57) 

(58) 


altimeter  bias 


altimeter  scale  factor 


ehsf 

c  =  static  pressure  coefficient 
sp 

=  barometric  time  delay 
6D  =  disturbance  input  to  vertical  channel. 

The  altimeter  bias  e  ,  more  commonly  known  as 
standard  setting  error  or  variation  in  height  of  a  constant 
pressure  surface,  varies  slowly  due  to  two  reasons. 

(1)  motion  of  the  vehicle  through  the  weather 
pattern; 

(2)  motion  of  the  weather  system. 

The  rms  variation  of  this  altitude  has  a  bounded  magnitude. 
In  this  thesis,  this  error  is  modeled  as  a  first  order  lag 
given  by 


e  =  -u  . .  e  +  w 
po  alt  po  po 


U) 


alt 


=  v/d 


alt 


N 


alt 


2ujalt  aalt 


where 


d  , ^  =  correlation  distance  of  the  weather 

alt 

°alt  = 


one-sigma  value  of  the  variation  of  altitude 
of  a  constant  pressure  surface 


V 


alt 


vehicle  speed 

power  spectral  density  of  the  white 

gaussian  noise  w 
r  po 


The  altimeter  scale  factor  e.  -  is  the  error  due  to 

hsf 

deviation  of  the  atmospheric  temperature  from  the  assumed 

temperature  profile  (Ref  7) .  The  indicated  altitude  error 

(e,.  )  is  of  the  form 

'temp 


'temp 


"  (ehsf}  (h) 


Thus  it  can  be  viewed  as  an  altimeter  scale  factor  error. 
This  error  varies  slowly  with  location  and  time,  and  it  is 
assumed  to  be  constant  over  a  typical  navigation  flight 
duration.  Thus  it  is  modeled  as  a  random  constant 


(63) 


with  appropriate  standard  deviation. 

The  altitude  indicated  by  the  barometer  is  based  on 

the  static  pressure.  The  latter  is  taken  from  the  pressure 

measurements  made  by  the  pitot-static  tube  in  the  vehicle. 

The  altimeter  error  e  due  to  the  erroneous  interference 

sp 

of  the  static  pressure  is  (Ref  5) 


(64) 
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rv 


The  coefficient  c  is  nearly  constant  with  altitude  and 

sp 

thus  it  is  modeled  again  as  a  random  constant 


with  appropriate  standard  deviation. 

The  barometric  time  delay  represents  the  time 

required  by  the  static  pressure  in  the  cavity  of  the  pres¬ 
sure  transducer  to  adjust  to  the  static  pressure  at  port 
by  flow  of  air  through  tubing  during  vehicle  maneuvers 
(Ref  6) .  This  time  constant  is  nearly  invariant  and  is 
modeled  as 


Tb  -  0  (66) 

with  a  given  standard  deviation. 

The  disturbance  input  <5D  is  the  cumulative  effect  of 
all  other  sources  of  error  which  may  influence  the  vertical 
channel  during  the  vehicle  climbs  and  descents.  This  process 
has  not  been  modeled  as  an  additional  error  state;  instead, 
it  is  treated  deterministically.  That  is,  it  is  fed  to  the 
vertical  channel  during  the  time  when  the  vehicle  performs 
a  descent.  The  magnitude  of  this  error  was  selected  as 
200  meters  based  on  discussions  with  the  sponsor  of  this 
thesis.  Thus 


66 


where  (t^  -  t^)  is  the  time  during  which  the  vehicle 
descends. 

Trajectory  Selection 

The  trajectory  generator  is  needed  to  give  position, 
velocity  and  specific  force  throughout  the  interval  of 
study.  General  trajectory  generation  programs  are  available; 
however,  in  this  thesis  it  was  decided  to  conserve  the 
computer  resources,  and  a  trajectory  of  a  great  circle  path 
was  generated  by  a  set  of  closed-form  expressions  for 
position,  velocity  and  specific  force  (Ref  12) . 

The  pattern  of  the  flight  path  selected  was  a  straight 
and  level  flight  at  600  miles  per  hour  at  11,000  feet  for 
a  duration  of  500  seconds,  followed  by  a  dive  at  the  rate 
of  6000  feet  per  minute  for  100  seconds  (with  a  corresponding 
decrease  in  ground  speed) ,  and  finally  leveling  at  1000 
feet.  The  altitude  profile  generated  is  shown  in  Figure  14, 
with  the  vertical  velocity  as  shown  in  Figure  15. 

The  mission  scenario  envisions  a  disturbance  input  to 
the  vertical  channel  at  t  =  500  seconds  just  at  the  time 
the  vehicle  starts  descending,  and  ends  at  t  =  600  seconds 
as  the  vehicle  levels  off.  At  time  t  =  610  seconds,  the 
vehicle  is  required  to  perform  a  TERCOM-update  ending  at 
time  t  =  660  seconds.  It  is  during  this  last  section  of 


Figure  14.  Altitude  Profile 


Time  (secs) 
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Figure  15.  Vertical  Velocity  Profile 


that  the  behavior  of  the  vertical  channel  is  analyzed.  The 
various  time  intervals  were  given  by  the  thesis  sponsor. 

Monte  Carlo  Simulation 

A  generalized  Monte  Carlo  simulation  program  (SO FE  - 
a  generalized  digital  simulation  for  optimal  filter  evalua¬ 
tion)  (Ref  13)  was  used  to  propagate  the  error  states  over 
the  total  time  interval  of  700  seconds  (covariance  analysis 
program  was  not  done  due  to  non-zero  mean  disturbance  input) . 
The  companion  plot  program,  SOFEPL,  was  used  for  generating 
subsequent  plots  (Ref  14) . 

SOFE  requires  both  the  true  error  states  and  the  filter 
error  states.  Since  there  were  no  filter  error  states  for 
this  work,  a  dummy  filter  error  state  was  programmed.  The 
various  user  input  routines  to  SOFE  for  one  set  of  vertical 
channel  gains  are  given  in  Appendix  C. 

The  integrator  in  the  basic  SOFE  cannot  handle  step 
changes.  Since  the  very  nature  of  this  thesis  involved 
step  changes  due  to  disturbance,  altitude  and  vertical 
velocity,  these  were  approximated  as  cosine  functions  over 
a  small  interval  of  time  (At  =  0.01  secs)  at  each  corner. 

This  device  enabled  the  integrator  to  work  properly  and 
made  no  substantial  impact  on  the  results. 

Three  plots  were  generated  for  each  flight.  Error 
states  3  (altitude  error) ,  6  (vertical  velocity  error)  and 
44  (barometric  error)  (see  Table  IV-1)  were  plotted  using 
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SOFEPL.  The  three  plots  generated  were  for  each  set  of 
vertical  loop  gains;  i.e.,  classical,  improved  and  combined 
gains.  Thirty  Monte  Carlo  runs  were  carried  out  for  each 
case  to  get  an  ensemble  average.  For  one  set  of  loop  gains, 
the  statistical  results  of  100  Monte  Carlo  runs  were  not 
significantly  different  from  the  results  of  thirty  Monte 
Carlo  runs.  To  conserve  computer  resources,  thirty  Monte 
Carlo  runs  were  used  in  all  the  following  analyses.  For 
further  insight,  one  Monte  Carlo  run  was  also  carried  out 
for  each  set  of  gains.  For  comparison  purposes,  the  starting 
point  of  the  random  number  generator  was  set  to  a  fixed 
value  for  all  three  sets  of  gains.  That  is,  each  separate 
set  was  generated  with  the  same  noise  realization.  The 
results  obtained  from  these  simulations  are  given  in  the 
following  chapter. 


V.  Results 

Basic  Cost  Function 

Recalling  from  Chapter  II,  the  basic  cost  function  was 
fc4 

J(K)  =  /  (5h  )  2  dt  (6 

fc3 

For  the  above  cost  function,  the  barometric  data  was  assumed 
perfect  with  no  uncertainties  whatsoever  (Figure  1) ,  except 
during  the  interval  (t2  -  t^)  (see  Figure  2)  at  which  time 
the  disturbance  was  fed  into  the  vertical  channel.  Equation 
(68)  was  minimized  through  the  search  routine,  and  the  optimum 
values  of  the  three  vertical  loop  gains  are  given  in  Table  V-l 

TABLE  V-l 

Optimized  Gains  (Basic  Cost  Function) 


Basic  Cost  Function  J (K)  =  /  (6h.)2  dt 

fc3 


Vertical 

Loop 

Gains 

Optimized  Result 

Units 

K1 

935678.67 

-l 

sec 

K2 

0.02 

-2 

sec 

K3 

0.0001 

_  3 

sec 
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The  message  from  Table  V-l  is  apparent.  With  the  cost 
function  of  Eq  (68) ,  the  vertical  loop  gain  needs  to 

be  set  to  essentially  infinity  regardless  of  the  values  of 
K2  and  •  If  we  analyze  the  set-up  of  the  cost  function 

in  more  detail,  the  result  is  not  surprising.  As  said 
earlier,  the  barometric  data  is  assumed  perfect  except  during 
the  disturbance  interval;  i.e.,  the  barometric  data  is 
perfect  before  and  after  the  disturbance.  The  optimum  INS 
altitude  estimate  after  the  disturbance  would  be  the  per¬ 
fect  baro  output,  and  tight  tracking  control  will  minimize 
the  INS  variation  from  this  ideal  baro  indication.  So,  the 
large  value  of  gain  can  be  anticipated  for  this  highly- 

simplified  model.  If  baro  altitude  were  not  perfect,  then 
the  gain  would  settle  for  a  far  lesser  value,  thereby 

indicating  that  the  inertial  system  does  not  truly  believe 
in  the  data  from  the  altimeter  due  to  uncertainties  in  the 
latter.  For  this  case,  however,  the  data  from  the  altimeter 
is  perfect,  thus  the  gain  must  be  set  to  infinity  to 

track  the  altimeter  without  any  lag.  A  constraint  optimiza¬ 
tion  routine  or  addition  of  a  term  in  the  cost  function  to 
penalize  huge  values  of  (or  K 2  or  )  could  be 

effectively  used  at  this  stage.  Such  a  procedure  is  incon¬ 
sistent  with  the  objectives  of  this  analysis  and  was  not 
pursued . 

The  optimized  values  of  the  three  gains  of  Table  V-l 
are  not  very  accurate  because  of  the  inherent  limitations 


of  the  integration  routine  used;  however,  these  values 
give  insight  into  the  behavior  of  the  vertical  channel. 

New  Cost  Function 

With  the  addition  of  uncertainties  in  the  vertical 
channel  and  barometric  data  (Fig.  4),  the  cost  function 
was 

fc4  - 

J (K)  =  /  (6h1)2  dt  +  (t4  -  t3)  [ (6h2) 2 ] 

fc3 

It  may  be  of  importance  to  note  that,  in  the  minimization 
of  Eq  (69),  a  weighting  factor,  3  ,  can  be  introduced 
such  that 

fc4  - 

J(K)  =  (6)  /  (6hL)2dt  +  (1-3) (t4-t3) [ (5h2) 2]  ' 

t3 

For  a  value  of  3  (between  0  and  1) ,  it  is  possible  to  have 
any  combination  of  the  mean  squared  error  due  to  the  dis¬ 
turbance  and  noise.  This  in  effect  scales  the  size  of  the 
deterministic  disturbance.  If  3  is  set  to  0.5,  then 
Eq  (70)  reverts  back  to  the  form  equivalent  to  Eq  (69) . 

To  provide  a  baseline  design  and  performance  against 
which  to  compare  the  optimized  performance,  the  classical 
set  of  gains  is 


(71) 


K2  =  3.0307  x  10-4  sec-2 

=  1.0  x  10~6  sec”3 

The  only  rationale  given  for  these  gains  is  that  the  third- 
order  control  system  has  a  triple  pole  with  a  100  second 
time  constant.  This  design  specification  allows  the  INS 
estimate  to  average  out  the  high  frequency  barometric  noise, 
but  is  not  optimal  in  any  sense.  It  may  also  be  interesting 
to  compare  the  performance  with  the  optimized  set  of  gains 
obtained  by  Widnall  and  Sinha  (Ref  4)  for  the  mean  squared 
velocity  error;  they  are 

=  1 . 004  sec" 1 

K2  =  4.17  x  10~3  sec-2 

K.j  =  4.39  x  10-6  sec-3 

As  stated  in  earlier  chapters,  the  magnitude  of  the 
disturbance  input  to  the  vertical  channel  was  assumed  to 
be  200  meters.  With  such  a  disturbance,  the  minimization 
of  Eq  (69)  led  to  the  optimized  set  of  gains  as  given  in 


(72) 


TABLE  V-2 

Optimized  Gains  (New  Cost  Function) 


New  Cost  Function 

fc4 

J (h)  =  /  (6hL)2dt  +  (t4-t3) [ (5h2) 2] 
t3 

Vertical 

Loop 

Gains 

Value 

Units 

K1 

0.631 

sec  1 

K2 

4.78  x  10"3 

-2 

sec 

K3 

6.335  x  10~5 

sec-3 

The  mean  squared  altitude  error,  with  the  values  of  the  noise 
spectral  densities  as  given  in  Chapter  II,  for  the  classical 
set  of  gains  (Eq  (71))  was  found  to  be 


Uh)*L 


=  595.879  m2 


=  (24.41) 2  m3 


The  corresponding  mean  squared  altitude  error  for  the 
improved  set  of  gains  of  Table  V-2  is 


In  calculating  the  mean  squared  error  as  given  in  Eqs  (73) 
and  (74) ,  the  mean  square  value  of  the  error  state  44 
(error  due  to  variation  in  altitude  of  a  constant  pressure 
surface)  (first  order  lag  shaping  filter.  Fig.  4)  was  sub¬ 
tracted  so  that  true  performance  comparison  could  be  made, 
between  the  classical  and  improved  gains.  The  mean  squared 
error  for  gains  of  Eq  (72)  was  not  calculated  for  reasons 
presented  later.  Accordingly,  (152.4  m) 2  was  subtracted 
(see  Chapter  II)  and  is  not  included  in  Eqs  (73)  and  (74) , 

This  performance  improvement  is  significant  relative 
to  the  classical  gains.  The  mean  squared  altitude  error 
is  70%  lower. 

The  poles  of  the  closed-loop  portion  of  the  vertical 
channel  are  the  three  roots  of  the  characteristic  equation 
(see  Chapter  II) 

S3  +  K^S2  +  (K2  -  2g/R) S  +  K3  =  0 

With  the  values  of  the  loop  gains  of  Table  V-2,  the  three 
poles  are  located  at 


P1  =  -0.6235  sec 

P2,P3  =  -3.75  x  10“3  ±  j  9.36  x  10-3  sec-1  (75) 

These  poles  have  a  time  constant  of 

=  1.604  sec 

x2,x2  ~  266.67  sec  (76) 

Comparing  with  the  classical  gains,  one  time  constant  is  a 
factor  of  100  faster?  the  other  two  time  constants  are  a 
factor  of  three  slower. 

The  individual  contributions  of  the  various  white 
noises  (see  Fig.  4)  to  the  mean  squared  altitude  error, 
for  the  noise  densities  as  given  in  Chapter  II,  are  shown 
in  Table  V-3. 

Table  V-3  shows  that  for  the  classical  gains,  the 
mean  squared  altitude  error  is  dominated  by  the  short 
correlation  time  acceleration  error  and  more  so  by  the 
altimeter  error  (first  order  lag) .  For  the  improved  gains, 
the  contribution  of  altimeter  error  and  short  correlation 
time  altimeter  error  is  the  greatest. 
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TABLE  V-3 


Contribution  of  White  Noises 
to  Mean  Squared  Altitude  Error 


Noise 

Density 

Mean  Squared  Altitude 

Error  (m) 2 

Classical  Gains 

Improved  Gains 

Qal 

15.0 

0.04064 

Qa2 

1.875 

1.686  x  10"3 

Qbl 

2.074 

31.937 

Qb2 

576.93 

22.82 

TOTAL 

595.879  =  (24. 41) 2 

54.899  =  (7 .409) 2 

Unfortunately,  under  the  presence  of  the  disturbance, 
comparison  cannot  be  made  between  the  optimized  gains  for 
the  mean  squared  velocity  error  found  by  Widnall  and  Sinha 
(Ref  4)  as  given  in  Eq  (72),  and  the  improved  gains  for  the 
mean  squared  altitude  error  of  Table  V-2.  In  addition,  the 
gains  of  Eq  (72)  are  the  optimized  gains  for  the  vertical 
velocity  error,  whereas  the  improved  gains  of  Table  V-2 
are  for  the  altitude  error  and  a  performance  comparison 
between  these  sets  of  gains  would  be  pointless.  To  gain 
insight  into  the  nature  of  the  optimized  gains,  the  disturb 
ance  was  set  to  zero,  and  using  the  values  of  the  dynamic 


9 


driving  noises  used  by  Widnall  and  Sinha  (Ref  4)  (with 


the  first  order  lag  for  noise  instead  of  random 

walk  as  done  in  Ref  4) ,  the  optimized  gains  obtained  are 
as  given  in  Table  V-4.  The  gains  of  Table  V-4  show  a 
considerable  departure  from  Table  V-2,  especially  the  gains 
and  .  In  essence,  as  pointed  out  by  Widnall  and 

Sinha  (Ref  4) ,  the  gain  primarily  depends  on  the 

strength  of  the  noise  sources  Qb2  and  Qbl  .  The  value 
of  gain  for  Case  II  is  in  excellent  agreement  with 

the  formula  (Ref  4) 


K 


1 


y  ioo 


sec 


0.51  sec 


The  gains  of  Table  V-4  (for  zero  disturbance)  are  not 
explicitly  required  for  this  thesis  since  the  very  objective 
of  the  thesis  was  to  find  optimum  gains  due  to  non-stochasti 
(disturbance)  and  stochastic  inputs.  To  obtain  further 
insight  into  the  nature  of  the  optimized  set  of  Table  V-2, 
it  was  necessary  to  find  their  sensitivity  to  the  time 
intervals  t^  ,  t2  ,  t^  ,  and  t^  (see  Fig.  2) ,  which 
is  analyzed  in  the  following  section. 

Sensitivity  Analysis 

It  may  be  recalled  from  Figure  2  and  Chapter  IV,  the 
vehicle  was  required  to  descent  for  a  period  of  100  seconds 
in  the  interval  (t2  -  t^)  ,  and  perform  a  TERCOM-update 
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TABLE  V-4 


<r* 


Optimized  Gains  for  Zero  Disturbance 


Optimized  Gains 

for  Disturbance  =  0 

Loop  Gains 

Optimized  Gains 

K1 

0.54 

sec  1 

K2 

7.77 

x  10-3  sec"2 

K3 

1.02 

,  -4  -3 

x  10  sec 

for  50  seconds  during  the  interval  (t4  -  t^)  .  Table  V-5 
shows  the  values  of  the  vertical  loop  gains  for  an  increase 
of  10%  in  each  of  the  time  intervals  (t2  -  t^)  ,  (t^  -  t2) 
and  (t^  -  t^)  ,  respectively.  Comparing  with  the  improved 
gains  (Table  V-2) ,  we  find  that  the  gains  are  very  sensitive 
to  the  time  interval  (t^  -  t2)  ;  i.e.,  after  the  descent 
and  before  the  TERCOM-update.  On  examining  it  more  closely, 
we  find  that,  increasing  the  time  interval  after  the  dis¬ 
turbance  interval  (t2  -  t^)  ,  it  is  logical  for  these  gains 
to  change  and  settle  on  the  steady-state  values,  because 
the  effect  of  the  disturbance  is  decreasing;  in  effect,  had 
not  the  gains  been  optimized  for  the  interval  (t^  -  t^)  , 
they  would  have  approached  the  values  as  given  in  Table  V-4. 


TABLE  V-5 


Sensitivity  of  Gains 


Sensitivity  Analysis 

Loop 

Gains 

Interval  (t^t^) 
increased  10% 

Interval  (t^-t^ 
increased  10% 

Interval  (t^-t^) 
increased  10% 

K1 

0.632  sec-1 

0.557  sec-1 

0.65  sec-1 

*2 

4.63x10  3  sec  2 

3.32xl0-3  sec-2 

4.9xlO-3  sec”2 

K3 

6.47xl0”5  sec-3 

5. 81x10” 5  sec”3 

5. 21x10” 5  sec”3 

It  is  natural  for  the  gains  to  approach  their  steady-state  value  in  the 
long  run,  once  the  effect  of  the  disturbance  is  over.  In  addition,  the 
very  slight  difference  between  the  intervals  -  t^)  ,  (t^  -  t^) 
and  those  of  Table  V-2  is  due  to  the  minimi  effect  of  the  disturbance. 

The  optimized  gains  for  the  disturbance  (Table  V-2)  were  checked 
out  in  the  simulated  flight  of  a  vehicle  performing  a  TERGQM-update 
and  the  results  are  shown  in  the  next  section. 

Simulation  Results 

As  stated  in  Chapter  IV,  one  flight  profile  was  used  with  three 
different  gain  sets;  classical,  improved  and  ocmbined.  In  the  first 
flight,  only  the  classical  gains  were  available  throughout  the  duration 
of  700  seconds,  and  in  the  second  flight,  the  improved  gains  of 


Table  V-2  were  programmed.  For  the  third  flight,  the  vertical 
channel  was  programmed  to  use  the  classical  gains  up  to  the 
time  t^  =  610  seconds,  and  then  switched  over  to  the  improved 
gains  of  Table  V-2  for  the  TERCOM-update  interval  of  (t^  -  t^) 
=  50  seconds,  and  finally  switched  back  to  the  classical 
gains  after  the  time  t^  (see  Fig.  2) .  The  gains  of  Table  V-2 
were  optimized  only  for  the  duration  of  the  TERCOM-update? 
therefore,  it  was  appropriate  to  program  them  only  for  this 
interval.  The  results  for  the  altitude  error  for  the  classi¬ 
cal,  improved  and  combined  flights  are  shown  in  Figures  16, 

17  and  18. 

On  examining  Figure  16  (classical  gains) ,  we  see  an 
initial  hump  around  time  t  =  100  seconds  with  the  error 
finally  settling  to  its  steady  state  value  at  time  t  =  400 
l  conds.  This  slow  rise  to  its  transient  peak  around  t  =  100 
seconds  is  due  to  the  inherent  lag  (r  =  100  secs)  in  the 
classical  gains.  Thus,  with  these  gains  it  takes  a  long  time 
to  build  up  the  error  and  settle  on  the  steady  state  value. 

At  time  t  =  500  seconds,  the  vehicle  went  into  a  dive,  and 
the  buildup  of  the  error  after  t  =  500  seconds  due  to  error 
in  altimeter  is  evident.  At  time  t  =  600  seconds,  the 
descent  of  the  vehicle  stops,  but  the  error  in  altitude 
takes  a  long  time  to  descend  down  and  follow  the  altimeter. 
Notice  at  time  t  =  610  seconds,  where  the  vehicle  «vas 
required  to  perform  the  TERCOM-update,  the  altitude  error 
is  still  in  excess  of  800  feet.  This  error  may  seem 


unrealistic,  but  for  the  kind  of  disturbance  and  the  flight 
trajectory  the  result  is  not  incorrect.  For  a  disturbance 
input  of  a  smaller  magnitude,  the  altitude  error  will  be 
correspondingly  lower  and  the  optimal  gains  would  be  differ¬ 
ent. 

For  the  case  of  Figure  17,  in  which  only  the  improved 
gains  were  programmed,  we  notice  a  considerable  change  from 
Figure  16.  Notice  how  quickly  the  error  builds  up  to  the 
steady  state  value;  this  fast  response  is  due  to  these 
gains.  The  curve  of  Figure  17  has  the  same  shape  of  Figure 
16,  except  that  it  has  sharp  response  features.  Also,  at 
time  t  =  610  seconds,  the  error  drops  sharply  since  the 
effect  of  the  disturbance  terminated  *-  time  t  =  600  seconds. 
The  altitude  error  at  the  time  t  =  610  seconds  is  less  than 
200  feet.  It  may  also  be  noticed  that,  while  in  Figure  17 
the  error  stays  at  a  constant  value  until  the  termination  of 
flight  at  t  =  700  seconds,  the  error  in  Figure  16  continues 
to  decrease  and  it  appears  to  go  to  zero.  This,  however,  is 
not  true.  Since  the  classical  gains  have  an  inherent  lag, 
the  system  is  going  through  a  transient  at  time  t  =  700 
seconds;  had  the  flight  duration  been  extended  beyond  700 
seconds,  the  error  would  have  risen  again. 

In  Figure  18,  we  see  a  combination  of  the  classical 
and  improved  gains.  The  altitude  error  follows  the  pattern 
of  Figure  16  (since  the  gains  are  the  same)  up  until  time 
t  =  610  seconds.  Here,  the  system  switches  over  to  the 


pattern  of  Figure  17,  and  maintains  a  steady  value  until 
time  t  =  660  seconds,  at  which  instant  the  classical  gains 
once  again  take  over  and  the  system  starts  following  the 
pattern  of  Figure  16  again. 

From  Figures  16  and  17,  at  time  t  =  610  seconds,  the 
error  drop  corresponds  to  an  improvement  of  about  70% 
which  is  the  same  as  stated  earlier  in  this  chapter. 

At  the  same  time,  the  plots  for  the  vertical  velocity 
error  for  the  three  flights  were  also  obtained  and  are  as 
shown  in  Figures  19,  20  and  21  for  the  classical,  improved 
and  combined  gains,  respectively.  As  before,  the  gains 
depict  a  similar  behavior.  For  the  classical,  the  error 
takes  a  long  time  to  build  up  and  decrease,  whereas  for  the 
improved  gains,  the  change  is  very  fast. 

The  plots  for  the  baromatric  error  for  the  three  cases 
is  also  shown  in  Figures  22,  23  and  24.  These  three  plots 
are  exactly  the  same,  thereby  confirming  the  fact  that  the 
random  number  generator  was  set  at  the  same  value  at  the 
start  of  these  flights  as  required. 

It  was  interesting  to  observe  the  behavior  of  the 
system  for  one  Monte  Carlo  case  out  of  thirty  runs. 

Figures  25  through  33  show  the  behavior  of  the  altitude, 
vertical  velocity  and  barometric  errors  for  the  classical, 
improved  and  combined  gains  for  one  Monte  Carlo  run,  all 
with  the  same  noise  realization. 


The  behavior  of  the  system  is  as  before,  except  for 
one  case.  Increased  noise  content  is  evident  in  the 
improved  gains.  Unfortunately,  this  behavior  is  typical; 
the  more  the  gain  is  increased,  the  more  noise  content 

appears  in  the  output.  Increasing  gain  allows  the  INS 

to  track  the  baro  altimeter  more  closely,  but  because  of 
the  noisy  contents  of  the  latter,  more  noise  is  apt  to 
appear  in  the  output.  However,  for  this  case,  during  the 
short  interval  of  the  TERCOM-update,  it  may  be  tolerable. 
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Figure  17.  Altitude  Error  (Improved) 
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.'igure  21.  vertical  Velocity  (combined) 
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Figure  23.  Barometric  Error  (Improved) 
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Figure  24.  Barometric  Error  (Combined) 
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Altitude  Error  -  One  Run  (Improved) 
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Figure  31.  Altitude  Error  -  One  Run  (Combined) 


Vertical  Velocity  Error  -  One  Run  (Combined) 


Figure  33.  Barometric  Error  -  One  Run  (Combined) 


VI.  Conclusions  and  Recommendations 


The  simulated  flight  has  demonstrated  that  using  the 
improved  gains  for  a  vehicle  carrying  out  a  TERCOM-update , 
lower  mean  squared  altitude  errors  are  possible  as  compared 
with  the  classical  gains.  Due  to  the  inherent  lag  in  the 
classical  gains  (time  constant  of  100  seconds) ,  these  gains 
become  unsuitable  for  such  a  mission.  Instead,  by  optimizing 
the  gains  for  the  period  of  TERCOM-update,  it  was  shown  in 
the  previous  chapter,  an  improvement  of  70%  was  achieved 
over  the  classical  gains.  The  classical  gains  with  their 
long  time  constant  have  an  advantage,  in  that  errors  build 
up  slowly;  however,  on  the  other  side,  the  errors  also 
decrease  slowly.  Thus,  the  performance  achievable  through 
the  classical  gains  is  at  its  best  for  a  level  flight  of 
long  duration.  Any  climbs  or  descents  of  the  vehicle 
degrade  the  performance  considerably.  The  optimized  gains 
with  a  fast  time  constant  showed  close  tracking  capabilities 
of  the  vertical  channel  to  the  altimeter. 

The  results  of  the  sensitivity  analysis  showed  that 
the  loop  gains  are  highly  dependent  on  the  time  intervals 
during  and  after  the  disturbance.  Any  change  of  more  than 
10%  on  the  time  intervals  (t2  -  t^)  and  (t4  -  t^)  would 
warrant  a  new  set  of  gains.  The  results  also  showed  the 
gains  to  be  highly  sensitive  to  the  time  interval  (t-j  -  t2)  . 
Thus,  if  this  thesis  needs  to  be  adapted  for  a  different 


set  of  time  intervals,  a  different  set  of  optimized  gains 
needs  to  be  searched  for  (using  the  search  routine).  In 
addition,  it  was  also  shown  that  the  gains  are  a  function 
of  the  magnitude  of  the  disturbance.  Consequently,  with 
a  different  magnitude  of  the  disturbance,  the  optimum  gains 
will  also  be  different. 

In  essence,  this  thesis  has  demonstrated  that,  for  a 
vehicle  carrying  out  a  TERCOM-update,  it  is  advantageous  to 
have  gains  which  are  radically  different  from  the  classical 
gains.  The  combination  of  the  improved  and  classical  gains 
showed  that  better  performance  is  achievable,  rather  than 
with  one  set  of  gains  only.  Although  the  optimized  gains 
show  greater  susceptibility  to  noise  than  the  classical 
gains,  the  effect  on  the  system  performance  is  not  disturbing 
due  to  the  short  time  operation  of  the  former  gains. 

For  the  adaptability  of  this  thesis  to  be  such  a 
requirement  in  the  real  world  environment,  it  is  necessary 
to  define  the  different  time  intervals  for  the  TERCOM- 
update,  and  then  to  calculate  the  gains  as  done  in  this 
paper.  The  combination  of  the  optimized  gains  with  the 
classical  gains  will  then  provide  a  lower  mean  squared 
altitude  error  and  enable  the  vehicle  to  carry  out  a  more 
accurate  TERCOM-update.  It  is  recommended  that  further 
study  be  carried  out  in  determining  the  optimal  combination 
of  the  classical  and  improved  set  of  gains.  The  time  at 
which  the  improved  set  of  gains  should  be  switched  into 


the  vertical  channel  needs  to  be  determined,  which  will 
give  optimum  balance  between  the  increased  noise  level 
content  and  optimal  performance  of  the  vertical  channel 
in  providing  lower  mean  squared  altitude  errors  at  TERCOM- 
update . 
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Appendix  A 

Instability  of  the  Vertical  Channel 


For  a  vertical  accelerometer  with  input  axis  along  the 
z-axis  in  the  local-level,  the  measured  specific  force  is 
given  by 


The  gravity  in  Eq  (A-l)  can  be  given  by  the  Taylor  series 
expansion  truncated  to  first  order 

g  =  go  "  T2  h  (A_2) 

e 

where 

gQ  =  gravity  at  distance  rQ 

R  =  radius  of  earth 

e 

h  =  height  above  earth. 
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Thus,  Eq  (A-l)  becomes 


6f 


2 


Sh  +  <$g 


or 

2g 

. « 

6h  -  -g—  <Sh  =  5f 

e 

Assuming  6f  to  be  constant,  the  solution  to  Eq  (A-4)  can 
be  written  as 


For  a  gravity  anomaly  of  10-6  g,  the  error  in  altitude 
after  one  hour  is  (Ref  15) 

<5h  =  28,900  ft 

This  shows  that  the  vertical  channel  is  unstable  for 
altitude  calculations,  and  external  altitude  information 
must  be  made  available  to  stabilize  it. 


(A-3) 


(A-4) 


(A-5) 


Appendix  B 


Minimization  Algorithm  Listinc 


I 


THE  OBJECT  OF  THIS  PROGRAM  IS  TO  MINIMIZE  A  COST  FUNCTION  WITH 
THREE  INPUT  VARIABLES. IT  UTILIZES  TWO  IMSL  ROUTINES  NAMELY  DCEAS 
AND  ZXM1V  TO  PERFOW  THE  JOB. THE  THREE  INPUT  VARIABLES  TO  THE 
Z!WIN  ROUTINE  ARE  THE  THREE  CAINS  OF  THE  VERTICLE  CHANNEL  OF  THE 
INS. THE  VERTICLE  CHANNEL  IS  MODELED  BY  A  SET  OF  FOUR  DIFFERENTIAL 
EQUATIONS  WHICH  ARE  SOLVED  PY  THE  DGEAR  ROUTINE. FOR  AN  INPUT  SET 
OF  CAINS  DCEAR  SOLVES  FOR  THE  COST  FUNCTION  AND  ROUTES  THE  RESULT 
BACK  TO  ZXMIN  WHICH  LOOPS  A  NEW  SET  OF  GAINS  .AND  THIS  PROCESS 
CONTINUES  TILL  CRITERIA  FOR  A  MINIM  UN  IS  MET. 


PROCRAM  TH8 

EXTERNAL  FUNCT 
CCMMON/DATA4/H(6  ) 
DIMENSION  C(3),W(R) 
INTEGER  M AXFN.N, 10PT 
REAL  K(3) 


INITIALIZE  INPUT  VARIABLES  FOR  ZXMIN. 

M-3 

KSIC-3 

MAXFN-500 

IOPT-3 

CAIN1-.82 

INITIAL  VALUE  FOR  VERTICAL  LOOP  CAINS. 

CAIN2-4.91E-3 

CAIN3-5.29E-5 

SET  SCALE  VALUES  FOR  LOOP  CAINS. 

K(1)-CAIN1/1.E-1 

K(2)-CAIN2*lOOO. 

K(3)-CAIN3*1 OOOOO. 

CALL  ZXMIN  TO  MINIMIZE  COST  FUNCTION. 


0  CALL  Z3W IN( FUNCT.M .NSIG.MAXFN, IOPT, K,H , G, F.W, IER ) 
IF(IER.EQ.O)CO  TO  20 
PRINTING  ERROR  MESSAGES. 


IF(  IER.EQ.  I29)THEN 

PRINT*, 'HESSIAN  NOT  POS.  DEF.  IER- '.IER 
END  IF 

IF (IER.EQ. I30)THEN 
PRINT*, 'IER-', IER 

PRINT*,  'MIN.  COULD  NOT  BE  ACHIEVED  TO  NSIG  DIGITS' 
REVERSE  SCALING 
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PR  1ST* ,  'CAIN1-'  ,F.(1  )*1  .E-l 
PRINT*, 'CAIN2-' ,K(2)/I000. 
pr  I  nt*  ,  ’  ca  i  n3-  • ,  k  { 3 )  / 1  ooooo . 

PRINT*, 'COST  FUNC-'.F 
END  IF 

IF(I£R.EQ.131)THEN 

PRINT* .  'MAXF'JN  EXCEEDED. . IER- • , IER 

CO  TO  10 

END  IF 

STOP 

PRINTING  MINIMIZED  VALUES  OF  VERTICAL  LOOP  CAINS 

0  PRINT*,'  ' 

PRINT*  ,'GAINl-',K(l)*l  .  F.-1 
PRINT*, 'CAIN2-' ,K(2)/1000. 

PRINT* , 'CAIN3- ' ,K( 3 ) /I 00000. 

PRINT*, 'COST  FUNC-'.F 
END 


THIS  SUBROUTINE  CALCULATES  COST  Fl'NC  AND  USES  DCEAR 
TO  SOLVE  THE  DIFFERENTIAL  EQUATIONS. 


SUBROUTINE  FUNCT(M,K,F) 

EXTERNAL  FCN.FCNJ 
DIMENSION  X(4),IUK(4),MK(60> 

.  CCMMCN/DATAl  /TIM EZ,TIME1,TIME2,,IME3,T 

IT  CCMM0N/DATA2  /DBARO,  A1 1 

CCMM0N/DATA3/CR  AV.RREC  ,CAIN1 ,  CAIN2 ,  GAIN3 
CCWMON/DATA4/H(6) 

COMM  ON /CE  A.R/ DIMM  Y  (  4  8  ) ,  SDIMM  Y  (  4  ) ,  I DWM  Y  ( 3  8  ) 

REAL  K(M),J1,J2 

SET  FIRST  ORDER  LAC  CONSTANT 
A-5.7931605E-4 

SET  SIOIA  VALUE  FOR  1ST  ORDER  MARKOV  PROCESS  (METERS) 
SICMA-152.4 

SET  NOISE  VALUES. 

QAl-2 . 4E-4 
QA2-1  .E-9 
QB1-100. 

QB2-2*A»(SI(*A**2) 

N-4 

SET  INITIAL  VALUES  FnR  INTEGRATION 
T*0.0 
X(l)-0 
X(2)-0 
X(3)-0 
X(4)-0 
T0L-2.E-10 
S-. 000001 
METU-2 
MITEP-2 
INDEX-l 
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TSTOF-160 


WRITE  VALUES  FOR  PLOTTING  ON  TAPE  5 
REWIND  5 
WRITE(5)  T, X 
TEND-O. 

ICOUNT-O 

CSAVITY  IN  METERS/ SEC**2 
CRAV-9. 80665 

RADIUS  OF  EARTH  IN  METERS. 
R-6378165.0 
RREC-l./R 


REVERSE  SCALING 

CAIN1-ABS(K(1))*1.E-1 

CAIN2-ABS(K(2))/1000. 

GAIN3-ABS(K(3))/lOOOOO. 

SET  TIME  VARIABLES 

TIMEZ-O. 

TIME1-I00. 

TIME2-110. 

TIME3-160. 

BETA-.  5 

SET  DISTURBANCE  M ACNITUDE (SQUARED )M ETERS 
ALFHA-2  00**2 


C-2  *CRAV*RREC 
CALL  UPDATE 

REINITIALIZE  DGEAR  FOR  STEP  INPUTS 
IF(T.EQ.100)THEN 
INDEX-1 
S-l.E-6 
END  IF 

IF(T.EQ.  1 10)THEN 
INDEX-1 
S-l.E-6 
END  IF 

IE(T.LT.TSTOP)THEN 

ICOUNT-ICOUNT+1 

TEND-T+5. 

CALL  OCEAR 
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CALL  DCEAR(N, FCN, FCNJ , T, S , X,TEND,TOL,METH, 

■H  ITER, INDEX, IWK.WK,  IER) 

PRINTISC  ERROR  MESSACES 

IF( IER. CT. 1 28 )THEN 
PR  I  NT*,  IER 
IF(  I ER • EQ . 1 32 )THEN 
DO  30  I-l.N 

CD-WK(IDIMNY(U)+I)/WK(I) 

PRINT*, CD 
CONTINUE 
END  IF 
STOP 

END  IF 

WRITER)  T.X 
GO  TO  5 
END  IF 

COMPUTE  COST  DUE  TO  DISTURBANCE 
J1-ALPHA*X(4) 

COMPUTE  VARIABLES  FOR  J2  COST  FUNC 
Bl-QAl/(2*((CAINl)*(CAIN2-C)-CAIN3)) 

B2“(CAIN1*QA2 ) /( 2  *  C (CA1N1*CAIN3)*(CAIN2-C)-(CAIN3**2) ) ) 
B3N-(CAIN1**2)*(GAIN2-C) 

B3N-B3N+(GAIN2**2 )-(CAINl*CAIN3) 
B3D-2*((CAIN2-C)*CAIN1-CAIN3) 

B3-B3N/P3D 

B3-B3*QB1 

B4N-(CAIN1*CAIN2)*((A*»2)*CAIN1+A*CAIN2+CAIN3) 

B4N“B4N-CAIN1*( (A**2 )*C*CAINl+(A**2 )*CAIN3+C*CAIN3) 

B4N-B4N+(A*GAIN2)**2-CAIN3**2 

B4D-A*GAIN1*CAIN2*((A**2)+A*CAIN1+CAIN2-2*C) 

B4D-B4D+A*C*CAINI*(C-(A**2)-A*GAINI) 

B4D-B4D-GAINI*CAIN3*(C+(A**2)) 

B4D“B4D+CAIN3*( A*C-A* A* A-CAIN3 ) 
B4D-B4D+GAIN2*GAIN3*(CAIN1-A) 

B4D-2*A*B4D 

B4-B4N/B4D 

B4«B4*QB2 


CCM PUTS  COST  DUE  TO  NOISE  INPUTS 
J2-(TIME3-TIME2)*(B1+B2+B3+B4) 

CQ1PUTE  TOTAL  COST 

F  -  (  BETAM1  >+ ( 1  *B  E  TA  )  *J  2 


PRINT  VARIABLES  ITERATIVELY 
PRINT*,  'GA1NT-'  .CAINT 
PRINT*,  'CAIN2-  *  ,GAIN2 
PRINT*, 'CAIN3-' ,CAIN3 
PRINT*, 'El-', B1 
PRINT*,  'B2-*  ,B2 

PRINT*, ’B3-',B3 
PRINT*, 'B4-*,B4 
PRINT*, 'Jl-’.Jl 
PRINT*, 'J2-',J2 
PRINT* ,  'COST  FUNC-',F 
PRINT*,'  * 

RETURN 

END 


THIS  SUBROUTINE  CALCULATES  SETS  THE  DISTURBANCE  INTERVAL 


SUBROUTINE  UPDATE 

CCMM0N/DATA1/T1ME2,TIME1,T1ME2,TIME3,T 

COHMON7DATA2/DFARO,A11 

IF(T.CE.TIME2  .AND.T.  LE.TWEI  )THEN 

DBARO-l 

All-O. 

ELSE  IF(T.CE.TI>!E2.AND.T.LE.TIME3)THEN 
DBARO-O. 

AU-l. 

ELSE 

DBARO-0. 

All-0. 

END  IF 
RETURN 
END 


THIS  SUBROUTINE  CALCULATES  DERIVATIVES  FOR  DCEAR 

SUBROUTINE  FCN'(N, T, X, XDOT) 

CO1M0N/DATA2/DB  ARO,  A1 1 
CCMM0N/DATA3/CRAV,RREC,CAINI,CAIN2,CAIN3 
DIMENSION  XCN).XDOT(N) 

XDOTfl )— X ( 2 )-((GAINl)*(X(l )-DBARO) ) 
XDOT(2)-((2*GRAV*RREC)-CAIN2)*X(l)+(CAIN2*DBAF.O)-X(3) 
XDOT(3 )-(GAIN3*X( 1 ) )-(CAlN3*DB ARO) 

XDOT(4 )-Al 1*(X(1 )**2 ) 

RETURN 

END 


THIS  SUBROUTINE  ACTS  AS  A  DtttMY  FOR  DCEAR 
SUBROUTINE  FCNJ(N,T,X,  PD) 

INTEGER  N 

REAL  X(N),PD(N,N),T 

RETURN 

END 
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Appendix  C 

SOFE  User  Input  Routines 


Introduction 


SCFE  (Ref  13)  is  a  Monte  Carlo  simulation  program  that 
helps  in  analyzing  integrated  systems  employing  Kalman 
Filter  estimation  techniques.  It  can  also  be  used  for 
propagating  the  navigation  error  equations  over  a  desired 
trajectory. 

SOFE  requires  both  truth  and  filter  model  state 
variables.  Since  this  thesis  had  only  truth  model  states, 
a  dummy  filter  state  was  introduced  to  satisfy  the  SOFE 
requirements.  Brief  discussions  on  each  of  the  input 
routines  used  in  this  thesis  are  described  in  the  following 
paragraphs.  For  additional  information,  Reference  13  can 
be  consulted. 


Tape  5 


This  file  is  the  input  to  SOFE  for  the  following 


information: 


(a)  Problem  title 

(b)  PRDATA  information  for  initializing  matrices 
and  time  intervals  in  basic  SOFE 

(c)  Initial  values  for  the  truth  model  error 
variables 

(d)  Input  for  the  USRIN  routine 

(e)  Plotting  information 


L8 


Tape  5  input  can  reside  on  card  decks  or  can  be  entered 
interactively . 

Subroutine  AMEND 

This  subroutine  is  used  to  apply  total  feedback  control 
after  a  certain  number  of  measurement  intervals.  Since  no 
filter  state  was  used,  this  routine  is  just  a  stub. 

Subroutine  ERDY 

This  subroutine  was  generated  to  calculate  the  non¬ 
zero  entries  of  the  9x9  fundamental  matrix  of  the  INS 
differential  equations.  This  subroutine  is  not  explicitly 
required  by  SOFE,  and  it  was  generated  for  the  purpose  of 
program  clarity. 

Subroutine  5STIX 

This  subroutine  is  required  for  computing  the  user 
defined  quantities.  Since  no  quantities  were  required  to 
be  computed,  this  routine  was  also  a  stub. 

Subroutine  FQGEN 

This  subroutine  is  required  for  the  filter  states 
which  were  not  present  in  this  thesis.  Thus,  this  routine 
was  also  a  stub. 

Subroutine  HRZ 

This  routine  is  again  required  for  the  filter  states 
and  thus  it  was  also  a  stub. 


Subroutine  NUUNIT 


This  routine  was  generated  to  convert  the  statistical 
input  to  computational  units.  Thus,  all  units  conversion 
calculations  were  done  here. 

Subroutine  SNOYS 

This  user  routine  adds  gaussian  random  samples  to 
specified  truth  states  to  simulate  the  accumulated  effect 
of  process  driving  noise  on  these  states  over  the  noise 
accumulation  interval. 

Subroutine  STABLE 

This  subroutine  was  generated  for  printing  a  table  of 
the  statistical  input  for  the  truth  model.  This  routine  is 
not  required  for  SOFE  explicitly?  however,  it  helps  in 
fault  analysis. 

Subroutine  TRAJ 

Since  no  external  trajectory  program  was  used,  this 
routine  was  generated  for  establishing  the  great  circle 
flight  trajectory. 

Subroutine  USRIN 

This  user  defined  routine  was  used  for  reading  and 


printing  input  data. 


Subroutine  XFDOT 


This  routine  is  required  for  the  filter  states,  and 
therefore  it  was  just  a  stub. 

Subroutine  XSDOT 

This  subroutine  contains  the  derivatives  of  the  truth 
model  and  it  also  initializes  the  various  error  states. 


100=LN1 5/VERT ICAL  CHANNEL  ERROR  PROPAGATION-CLASSICAL  VS  IMPROVED 
110=  IPRDAT  A 

120=  NF  =  1 , NS=50,N=0,NZF=0,NXTJ=?, 

130=  LXTJ=.F. , 

H0=  T0=0. ,  TF=700.f 

150=  DTPRNT=5.,  DTCCPL*5., 

UO=  DTN0TS=3. , 

170=  DTPRPL=10. , 

180=  LPRXF=  .F . ,  LPRDG=.F., 

190=  IPRRUN=1 ,  IPGSIZ=55, 

200=  LPP=.T.,  LCC=.T. , 

210=  T0LER=.O0O1  ,  HMAX=60.,HNIN=.0001, 

220=  ISEEP=~2361 268,  IPASS=30, 

230=  I 
240=  50*0. 

250=0. 

260=1 ,1,1. 

270=0,0,0. 

280=  ISIGOS  BYNAHIC( 1 )*10*0. , 

290=  GB (1  )=2*. 0005, .0007, 6F  < 1  )=6*.0005,GFF(I )'3*0. , 

300=  GSF(1)=3*5.,GH(1 ) =6*2.5, 

310=  AB ( 1 ) =3*25. , ASF <1)=3*25.,AN( 11=6*5.15, 

320=  BAR0=50. ,GDE=26. ,GND= 17. ,8A=35., 

330=  BASF=.03,SPC=1 .54E-04,AL=.25, 

340  =  1 

350=  ISTATS  DALTS=250.,D6RAVE=10.,DGRAVN=10., 

360=  DGA=60.,GYNDS(1)=3*Q.0,ACHDS(l >=3*0.0, 

370=  BASIGS=250. ,EDSIGS=26. ,DNSI65=17. ,GASI6S=35., 

380=  • 

390=  ICONTRL  LFDBK=.F. ,LINTOO= .F . ,KXSGO=1 , 

400=  i 
410=  1.0 
420=  3,0, 0,1. 

430=  6,0,0, 1. 

440=44,0,0,1. 

450=  0,0, 0,0. 

460=  TIHE (SEC ) 

470=ALTITUD£  ERROR 
480=P0SIT ION  ERROR  •FEET* 

490=V£RTICAl  VELOCITY  ERROR 
500=VEL0CITT  ERROR  *FPS* 

51 0=BAR0  ALTIMETER  ERROR 
520=ALTITUDE  *FEET* 


SUBROUTINE  A«ENO 


nn  4  optpi 


FTN  (.J.5S* 


11/08/32  15.33 


1 


5 


19 


15 


23 


25 


33 


SUBR  I'JTINE  A“£N0<IRUN.T,NP.XS.VXTJ,XF,XS.XTRAJ» 

c 

c  user-vritten  subroutine  ro  apply  total  feeoback  control 
c 

CI-XCN  /C%TRL/LF08K,LINT03.<XSSO 
0!  SESSION  XFINFt  . XSINSl.XTRA  Jf  SXTjt 

c 

LCGI  cal  LFD8K 

c 

return 

IF  {.NOT.  LF08K1  return 

c 

09  ICO  I  Pi.  3 

X  S ( I >  p  «S(I »  -  IF(I) 

«Ft ; i  p  a . 

100  CONTINUE 

03  2 C 0  IP  11.13 

XSII)  p  «S(I*  -  XFIIl 
xF< I )  p  9  . 

2C0  CONTINUE 

X  5  (  3  3  I  =  XS<33>  -  XF113* 

XS<«5»  P  XS<»>d  -  XF<1*> 

X  S  (  3  ?  I  p  XS(3?»  -  XFI151 
XFdCl  p  0. 

XF  (  1  3  I  ;  0. 

XFU51  P  0. 

RETURN 

C 

c  . 

ENTRY  A«ENOC 

RETURN 

£-.0 


CAPO  NR.  SEVERITY  OETAILS  OIAGV03IS  OF  P»OBLE« 

n  t  there  is  no  path  to  this  statement. 
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subroutine  e*ot 


74/74  C?7rl 


ftn  4.»*«44 


IlZ'S/ii  35.33 


33 

S 

i 

1 


W 


t 


1  7 


13 


SU34.UTI  M  E807INTAAJ.T7AJ.r?733 
C 

C  £40  7  CC“?UT£S  NON-EE70  ENT3I£S  F0»  7  HC  4  73  £  UB9L  OCK  THAT 

c  **«£$  up  the  fu.oancntal  navigation  errc*  ots»«:cs 
c  £/7t  also  pills  cr«niN  area  *sha*e*  with  :us8ent  values  of 
c  several  cop-'utco  vapjasles.  inclosing  i-t-e  velocj’ics. 

C  SPECIFIC  PC  8  CE  £  **.0  CCH«A\DED  AN  G  J L  A 4  VELOCITIES.  PLUS 
C  V  AN  0£  4  •  ANGL  £  VALUES.  'SHAPE*  VAtliSLES  »•£  USED  IN  iSDOT. 

C  7/037  A-.O  F5GEN  PI8  *  VA«I£TT  IP  StATE  ANO  COVARIANCE 
c  ce«auTAf res:. 

c 

CjPNCN  /EA«*-/G“£GA,SEl.E3a.G££ 

cj»»?n  /PAi-za-8 p: . -alp = i .»! • two®: ,  ®®o 

Cj»*:n  /V0A«P/C*1 .C*2.C«3 

c:  ■•■in  /:  ha  u/ vt, vt.  ve. vg, 87, ft. pe.wc  a, vct.wce.alfa.calpa.'alf  a 
OIpE-.SIGn  I--aj<n7*aj>,f».»<41> 


c  »t:«  cur 

QU!»CO 

3  L  AT 

=  •  ®ij( 1 1 

23 

V7 

=  *AJ(2> 

V  7 

s  *  *  A  J  C  J  1 

¥2 

= 

P« 

=  ’CAJfSl 

P  7 

s 

25 

ft 

=  :uj(7i 

ALFA 

2  rOAJI*> 

VAAIA9LC: 


C 

C 

C3»»utc 

8  SnA-  CO*  Aso  7t»P0»A»t  WA3  !  ABLES 

CALPA 

3  cssmlfai 

33 

SALFA 

*  :in(alpa» 

VC 

3  CALPA. V7  -  SALPA.V7 

VN 

3  54LP4.V7  .  CALPA. V7 

VG 

s  53PT(V»..2  .  V7. . 2 1 

PE 

3  CALPA. P7  -  SALPA.PT 

35 

PN 

3  SALPA.P*  .  calfa.pt 

CLAT 

3  c:s«  =  lat> 

Slat 

3  SIN«‘LAT» 

tl»t 

3  SLAT  /  CLAT 

3cla* 

3  1  . /CL  AT 

43 

;«cgn 

3  C“CGA.CLAT 

C«CG2 

3  C."EGA«SlAT 

‘HCt 

S-VN.P ‘EQ 

3hON 

3  VE.P'EQ 

4H0E 

3  VE.P'EG'TLAT 

45 

WE 

3  -  hOE 

UN 

s  8hCN.C«EGN 

WE 

3  8  hOE  *C  ®£  G  2 

WCI 

3  calfa. phgc  »  salfa.un 

WC  7 

3-SALFA.4hCE  ♦  CALPA. UN 

50 

WC  z 

3  CCG2 

«  42 

3  VE*4'E8 

C 

c 

EVALUATE 

the  t:»c-oe®eloe«.t  v:n-ee»o  eleven's 

c 

IN  tH[  PUNOAUEN'AL  C»®C*  D7NA»:cS  «4T»I7 

33 

iii  3  3“oe.°cl*t 

<21  3  2..«:»EGN.VN»C,,EGE*VE)*ah'N»V..4CLAT.4CLA’ 

(  37  S-2  ..;«CSN.VE-»-l3N.VE.4CLAT.»CLAT 

# 


r** 

k\  ‘ 

w  -  *  _  '  • 

- -  ..  . 

k 

hV 

r>„ 

[:• 

r.- 

- 

i 

■ 

.  - 

;: 

su3sour: 

»,£  £c  Of  7«/7«  OPfM  FTN 

1 1/09/32  15.55. 

v_ 

F9X9  (4) 

=-2.*:*EG2»¥C 

Fairs  (5) 

x-C“E62 

. 

63 

F9x9  c  s 

=  y  N*  Rh  0? • T  L  AT 

F9X9  (  7  > 

:-aH0;,«S^E3*»CL»T 

-  ' 

F9X9  C  9  ) 

:  ShOE<T"E3 

- 

F  9  <  9  ( 9  ) 

:-CXl 

F9X9U0  ) 

:  kh02**H0E*»hC5* X< 2 

-  * 

55 

F9I9C11 > 

=  shCN«s h02-’h;E» XX 2 

• 

F 91*  < l 2  ) 

r  2.«GCE*XFE3-<SH3V»XH3>(»7M0E*>H0C)-CX2 

1 

F  ?  X  9 (13) 

=-chOE*RPE0 

F?  X9 ( 14 ) 

=-rhcn»rreo 

F9X9 (15) 

:.iW02*a,E3 

:  - 

73 

F9X9< 1 6) 

=  S  =  CQ.RCL  AT 

FaX«<171 

i.jhCE*  TUIT-XK2 

f5x$U3» 

-•i  .<u2 

F  9  X9  (19) 

=  2.«WS 

F9X9(21) 

:  TL»T.a*E0 

75 

F9*9 (2  3) 

=  m2*C“ES2 

* 

F*X9 (29) 

=  -x«  i 

F9x9 (25) 

:-2  •  •-H3E 

F  9  X  9  (  2  9  ) 

=  -(WN»:,«E6N) 

F9X9 (29) 

=  c  hCE 

n 

F9xS  (  JC  » 

=  F  2 

F«x5<JU 

s-FN 

F9XJ<32> 

X-K2 

F? X?  ( 331 

s 

.- 

F5X5  <  J«l 

=  -F2 

arm 

55 

F?XS«J5> 

s  FE 

F? X 9  (  36) 

-  «2 

F9X*»(  57) 

5-ME 

- 

F9X«  <3»l 

r  F\ 

:- 

F9x9<  39) 

=  -FE 

=  3 

F9X3 (40  > 

-  -  W  N 

F9x7(«l) 

*  xE 

acru-N 

i  _ 

c 

c  . . 

• 

?5 

ENTRY  EROTO 

c 

C  EV«LU*TC  THE 

r:«C-:N0EPES3E*(r  VCV-2C70  £(.£**£  VTS 

C  15  the  FU*.0*«EST»l.  EX®CX  0T5»«ICS  »»T»Ix 

k 

3  REQ 

s  i./fica 

\ 

135 

F9 X9 (20  ) 

=  -£Q 

", 

F9x9(2  2) 

i  x’Ea 

‘ 

F7X9(26) 

:-«sEl 

F9X9 ( 2  7) 

=  1. 

-4 

SETUPS 

3  35 

C-iO 

4 

1 

— . 

I 

125 

.  * 

1 

L, 

.  „ 

-J i 

SUBOUTINC  ESn* 


7*/ 7* 


3T=  l 


FTN  *.>1»5S* 


11/03/12  15.33 


► 


1 


SUBROUTINE  eSTIX<:«U^*TtKP*'<S*'<<Tj*KF*XS.XT«4J*hT®*Pf> 

£NT*T  ESTIXO 
3  etu3* 

Z\Q 
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I'd 


SU993JTIKE  FOSE* 


7*/7% 


0O-2J 


F?V  «.P*56« 


12/C9/J2  15.33 


l  SUBflCUrisE  FQGESC  rPUSt  '  tSF.^SfV  *T  J*KF,  KSt  «T»*  F*OF  > 

C 

Q l “ENSIGN  XF<NF»**$<NS>.lT<»*J<N*TJ>,F(N2F>tQF<N?'J> 

■iETua** 

5  ENT*T  FOGESO 

CETU°N 
E‘«  0 


127 


3UB«:urist 


74/74  33  T=  1 


FT«|  4,!*H4 


11/34/12  15.SJ.J 


iUB»CuriS£  «»2C:»U1.r ,NF,4t.S*TJ,E0*,TO«.T*f NTc.PF, 


»£TU3* 
EN747  H8/0 
FE7U-X 
EMJ 
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3UBS3UTl*iE  NUUNIT  JW?«  0PT:1  FTS  *.B»5S*  11/05/32  15.33 


t  subssutine  nuum' 

c 

C  NlW'lIt  CONVERTS  STATISTICAL  IN®UT  '3  C3«PUTATICNAL  UNITS 

c 

5  C2*“CN  /CJS®S/0ALTS,0G®AVE.3G®AV';,D6A 

C3«**CN  /»Ari/Q»®P!."ALFPI  .®!.H0®l.®PO 
CC "NC  N  /NOIS3S/STN3St3».ACN3'«3» 

C3““3N  /SIGO/OTNA«ICtl3».SSC3»,GF  <S1.GFF<3>.GSF<3>.6«(6I  t 
•  ABC  3  I  .  CF  <3  I  .  »■<  5)  .  BA  CO  .  GDE  .GNO.GA, 

13  •  BASF.3PC.AL 

CO ■NON  /I  I G»A S/3 A  SI  63. COS IG 3. 3  NS  I G3.6AS1 GS 

c 

c 

C  FC®<«  C2NVE3SI0N  FACTi"S  F  0*  CHANGING  INPUT  0  AT  A 

15  c  to  co*putat:;nal  units 

c 

C  •••  »PO  :  OAOIAnS  PE®  DEGREE 

c  •••  SPH  =  SECCNOS  PE®  «Cu® 

C  •••  ACCPG  :  FT/SCC/3CC  PE®  SEE 

21  C  •••  ACCFU6  :  FT/SEC/SEC  PE®  “IC®3  SEE 

C  •••  SPAS  :  ® A3! AnS  PE®  A ® C  5ECCND 

C  ...  FPN«  i  FEET  ®E®  NAUTICAL  »!.E 

C 

=PO  :  mAlF®I/®0. 

25  :ph  :  3 6C  C  « 

ACCPG  :  32.2 

ACCPLG  :  32.2/IOOOCC3. 

®PAS  :  ®P0/35C3. 

FPNH  :  1552  ./.  33®! 

31  C 

C  change  STATISTICAL  INPUT  DATA  '0  1 3 '•PU-T  A  T I  3N  AL  UNI'S 
C 

03  ic  :=r,9 

1C  OTNA«IC«:»  :  QTNA-ICO/133’. 

35  DC  23  1:1.3 

GB(U  :  GB(I»*PPO/S»H 

ABC  I)  s  ABC  )*ACC®UG 

GTNOSCI)  =  <«GTS0S«!I*®P3/Ss,O»*2.)/S®N 

GFF  ( !  I  :  GFFU  >•  »PD/ <  S»H.  ACCPG*  AC  CFO  > 

*3  GIF  (1 1  :  GSFCI/13CC  33  3  . 

A  SF  1 1  I  :  AIFII »/10::i33. 

acnos  cii  =  <iacnos«h»ac:pjj»*»2.»/sph 
2C  CONTINUE 

c 

*3  OS  30  1:1.6 

GFTII  :  GF«t».»PO/«SP“»*CCPS» 

G»(I1  :  G««tf«PA5 
A“<  I  »  :  A«( I  I • ® P  AS 
30  CONTINUE 

53  C 

OAtTS  :  0ALtS*FPNN 

0  G®  A  VE  :  OGPA»E*FPNP» 

OG®  A  VN  :  3G® AVN.FPN* 

OGA  :  OGA.FPNN 

55  GOC  :  GOE  *  ACC  PuG 

GNO  :  GVO.ACCPUG 

GA  :  GA.ACCPUG 
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su03out:?.e  snots 


74/7%  CP  T  s  1 


FTS  4.?«5S4 


ll/C«/32  15.53 


nr 


t 


5 


11 


15 


21 


25 


11 


35 


O 


C 

C 

c 

c 

c 

c 


e 


c 


c 

c 


subroutine  s‘.ovs<irun>t>n;>vs>v»t j>  »f . xs .  *  rc  a j» 


USER-VPITTEN  SUHAOUT  1NE. 

*033  GAUSSIAN  =AS00N  SAMPLES  TO  SPECIFIED  T  =  UTM  states  '0 
SIMULATE  fit  ACCUMULATED  EFFECT  IF  PROCESS  0»!V!NG  f.OISE 
OS  THESE  STATES  OVER  The  NOISE  ACCUMULATION  INTERVAL  OT. 


COMMON  /C0K5S/0ALTS.0GRAVE.3GR  AVN.OSA 
C3"MCN  /S0IE3S/C»SDS(3».ACNDE«J1 
COMMON  /SIOMAS/8ASIGS.EOCIS3.DNSISS.SASISS 
OIMESSION  lF<  NFl t XS<5S>  t ITRA J( N«T Jt 


OT  s  T-'OLO 

VE  =  XT»AU<21 

VN  S  *TP*J<J( 

VG  p  SQ»T(ve*VE*VN*V*,> 


SIGH 

«S(lll 

SIG12 

«S<12) 

CIG1  3 

XST131 

SI  S3  2 

XSI32I 

SIGH 

»S(33» 

0163a 

»SU«1 

SIGMA 

*s<»*> 

:igm5 

«S(«5I 
SIG*S 
XSIMSI 
:IGM» 
XS (A  7 1 

tolo 

PET  U" N 


S0PT«0T«GTSDS(1»» 

XSU1  CSAUSS<O..SI  HI  » 

SORT!  0T>6T'.0SC2)T 
XS(12 I •GAUSS* 5.. S  IG121 
SOA  T( 0T*6TN0S( 3>> 

AS<13t*GAUSS(l. >315131 
SORT<OT.ACSOS(1  >» 

»S< 32 >»G*USS(C . .S IS32 > 
G0»T(0T*AC“.0SI2I» 

XS«33»*GAUSS(0. .515331 

sor  r«  ot  •  ACf.os  c ;  n 

<s(5a)«gauss(:.ii:s3a> 

8AS;GS>S0PTtl>-EAPt»?.*3T*VG/3ALrSl> 

xS<«A)«GAU5SO.t5tG4»> 

eos  :  si  -sort  a  .-e*pi-2..ot>vg/ogpave»* 

*S«A5».GAUS3«1.>SI JA5> 
OnS:gS«3QPT<1.-Exp«-2.>OT.vgtog5»vni » 
XS<A6I»GAUSS<0«>5 I  Gas  ) 
GA3;GS»SaPMi.-EXPC-2.OT>VG/0GA  )) 

<S( AT)«6AUSS(0.>S1SA7) 

T 


Entrt  ssorso 
tolo  =  T 

RETURN 

CM) 


Aj 


3U830UTI  VC  STABLE 


7A/74  0°  T  =  1 


ft* 


11/3S/12  15.33 


3ub»:ut:ne  stable 

stable  faints  a  tasle  of  stati sn :al  in»ut  oat*  fop  the  t»lith  model 

CO-XCN  /COA»S/DALfS ,06°  AVE ,00  A  AWN, OCA 
CiAHCN  /£4»’N/0*£GA,aEG,ES8«GEE 

cop«*-:n  /voi oos/stnosi 3> .aonos i « » 

cjkacn  /sigpas/sasigs.eosi ss.onsigs.gasiss 

CO“«CN  /S  IGO/OTna-ICUO  I  .681*  I.GFlii  ,GFF  (  3  *  .  GSF  <  J  >  ,  6«<  *  I  . 

•  A8( ST  *  ASF (31 ,  A*  ( 3 »  ,  BAXQ ,  GDE «  GNO  « 6A , 

•  BASF.SPC.AL 

01  *£N3I3N  o:1  ECT  t  3  >  , 30U»CEt l*  )  ,  AA  IS  <  1 2  '  •  Ax«!  S <61  ,OTNA*(  13  * 

DATA  AXIS/**  *,*X  *,*T  *.*T  *.*2  *.*2  *«*  INFUTTX)*,*  SPINTTl*, 

•  •  sp:n<x»*,*  if'®uTtr>»,»  s»is(r»»,«  :n»ut<2>*/ 

DATA  AX»:S/*ABT  T*,*ABT  7", •AST  X*.*A8T  E*,*ABT  X*,*ABT  Y*/ 

OATA  OHECT/**  *.*Y  *,*2  •/ 

DATA  0YNAN/*,CC3TLTIFT*,*  FEET*,*  FEET*,*  FT/EEC*. *  FTF3EC*. 

•*  FT/3EC*,*  AO* ,  •  "PAD*,*  «»AD*,*  FT/3EC**2*F 

OATA  S3U-CE/*EAST  *, *LCNG!TjDE*,*NO»TM*, "LATITUDE*. *ALTI ruDE*,* 

•  ,*ea:t*,*wel  :cirr*,*>.GAT^*,*vELc:r  t  r*.*vEaT:cAL*.*vELCci  *t*. 

•  •ea5t*,*att:tude*,*'<ortm*,*att:tuoe*,*wea»!cac*,*atfituoe*/ 

851  TE  table  OF  TAL'Th  MCOEL  SYATIS'IDAL  INPUT  OA*A 
VMITEIG.IOO  I 


W*ITE<6.0C0  1  1  ,SCU»C£«1 ».S3J"CE«2 1  .DynA"! C < l > *AE 0. DYNA«T 1 » 
MAI TE 16.2  00  I  2 ,3CU»C£« J) .SOJACEC*  > .DYNAM! C C 2  * *R£ 0.07 NA«1 2 > 
t!  =5 

00  13  :  =  3,b 

8"  I  TE  IS, MS  >  :  .SCUBCECI  >  .S3yACS«N*,0TNA«:CTt  »  ,0»‘.A»T  :» 
r:  =  i:  *2 
continue 


JAITET6.0Q0T  OYNAPICIIO  »,3TNA«U3> 

03  20  1x11,13 
LxI-13 

8»I  te (4,*:o )  I f0i3ECT<N>,6B(N>,GYN0S(S> 
CONTINUE 

00  30  I  x  1  A ,  1  f 
N=:-i3 
N'.  =  N»6 

8s  I  T£  (  6 , 3C  0  >  1,Ax:s(N>>axIS(nn»,GF(N) 
CONT : NUE 

00  «0  1:20,22 

\xt-l-* 

u«trc(i,t:oi  :  ,oi»ectin»,gffi  ;» 
continue 


0:  50  t : 2  3 ,2  3 
N  =  I  -  2  2 

«=ITE(i,7.C>  1 ,OI’ECT«N»,GSF(n» 
CONT : NUE 


'!* 


su8»aurisc  stable 


7*/T*  OPI-l 


ft s  *.!*«s* 


11/11712  15 


nr 


63 


65 


71 


75 


53 


65 


•1 


75 


103 


105 


113 


C 


60 

C 


70 

C 


50 

C 


1C 

c 


c 

103 


3C  3 
3C0 

«CO 

5C3 
6  C  0 
703 

■ICO 

703 

1030 

HOC 

1333 


O:  60  1=36.31 
=1-25 

Wc  I  TC  l  6  C  0  )  l.*«IS16>.»«*I31'l>.GH4» 

CONTINUE 

01  70  1=32.3* 

5=1-31 

u«1tc(6.ico>  i  .di’ECtii) ,*b is > .ac  .osini 
cost:  sue 

00  §3  1=35.37 
5=1-3* 

w-irc(6.i30o>  iiOurriM.i  Fm 

CONTINUE 

00  13  1=33. *3 
5=1-37 

«»;te(6.iioo»  i.»*is<5».«»“t>(6i.*«is» 

COST  0  5UE 

9*505  =  2.* HASiGS.BAS; 63/0*175 

Ei)  1505  =  2.»E05:GS‘EDSIG5/00F**E 
OnnOI  =  2..055ISS*05$16S/OG»*1,t 
6*505  =  2«*GASlGS*GAStGC/3SA 

*3  I T E (  6 , 1200  7  BAFO.BANOS.IAL'S.GIE.EO^NOS.OG^IFC.GND.ONNOS. 
DG»»»5,GA.G*10S.06A,B*SF,;»C,»L 


aETU55 

F0*«*r«*l*//-56.*T*UE  E»*?1  5 7 *7  I  671 CS*/ /T65 .2*1 • .• « . 7= J  . 

•  *:5Pur  P»=*«ETE»S  •,22<*.*>//T5,*T1UT..4,T->0,*50ISE*/ 

•  T5,*S*ATE*.T‘5.*IM't*L  1 5  £■  ,  r  7 1  ,  *SPECT*  »L  •  .71 1  7.  *C  07  7  EL*’  I  C5* , 

.  /7S,*:soei*.ri j.*ei*c*  s'u«:r*,r«*.*Eiic»  pcoel*. 

•  Tl5.*tIG*»  V*LUE*»T70,*DE5SI7t*.T11T,*p*;»*hE’ER*/1 
F0R1AT177, 12, 713, 2*7.71*. *D76A»;C*. *62, GI2.1.A10I 
Foa“»T(77,*u*.ri;,*vE7Trc»L  »:cELE7»nc6»,r**.*0T5*«ic». 

•  762 .612. *.*101 

F0m*7<77,12.7li,»2,T15.»6T»a  01 IFT** 7*4. •» » 5005  WALK*. 

>  T62.612.*.*  OEG/h4*.T»1,*C*,610.2. 

•  •  OEG/mR  |. *2/51*1 

F0R1AT1T  7,I2,713,A2,*GT*C  G-SC1S  01 1 F7  ,*,  *1. 7*«,  •«*500»  C0NSTA5**, 

•  762,612.*.*  0E6/H* /G*  7 

F 0***7177, 12, 713, *2.*G7*0  G*6-5ElS  0»IF7  *.724,7**, 

•  •=  *50 C 5  CONST  *51 *.r62, 312*4.*  0EG/H» / 1 G*G I • » 

F3*"»7|T7, 12.7U. *2, 715.*ST10  SCALE  F*C7C1* , 74* • *NOC»  C05574N'*, 

>  7  =  2, G12.*.*  PP4*7 

FQ1»*7«77,12,713.*2.*Gr»0  H  S  *L  I G  N*»C5  7  *,  *6. 7*  »  ,  *4*5005  CONSTAN**, 
7*2.012. ♦. *  A’C  EEC* • 

FC«i»7«7 7,12, 713«»2. 715. •*CCELE,3“E7E1  81 *S* »T»*,***NOC1  WALK*. 

I62.G12,*.*  UGCE*.7=?,*(*.SI3.2.*  gGEE >• • 2 /h* • | 

FS»»*717 r,12,7l3.*2.7l«,**CCELE40»E7C7  SCALE  FACTOR*. T*4, 

•  =  »50:*  C  )5ST  467  *. 762.612  . »,*  P  =  **7 

FO»»*l «7 7,12,71 J, *2, **CCEL  *1 5*L l G5«E1T  *,*6,744, 

•  =  *50  C5  c:5j7*57*.762,GI2.4,*  *®C  SEC*7 
FC°i*7 (7  7 ,**•• ,71 3»*8*®0  *l7 I *E  7  £  1  8;*S*,t*4, 

•F:»;r  c«ot»  «»**:»*. *s2.Gi2.*,*  fee7*,tig.gio.2, 

•  FEE'..;/5i*,7lt2,Gl0.3.*  1*UT  *ile:*/ 
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sua»3ur:sE  stable 


m/m  c»r  =  i 


ftn  ,.g.5s* 


ll/CSM 


5 


23 


23 


•  r;,*«*»,Tij,»E  g»av:tt  3Eflecti3*,‘,m«, 

•  ‘FIRST  C'DER  «AR*CV‘,TS2,G12.».‘  UGCE ‘ . T5 C ,61 3 .2, 

•  •  UG££*,2/*i“.ni2,G13.2,‘  S  A  U  T  ‘ILES*/ 

•  r».*«6*.m.*N  GSAWITT  3EFLECT  ICN‘,M*f 

•  ‘FIRST  3»DE»  «ARACV‘,TS2,S1>.*,‘  uGEE‘ ,T«C ,61 8 .2 , 

•  •  uGEE,,2/s***.rii2,ci:.2.*s»uT  «:les’/ 

•  T7.‘A7‘,TlJ.‘GSAV:Tr  ALS-LT* . T»A , 

•  ‘FIRST  Ca  0£R  FAA*CV*,’62,S12.*,*  UGIt^fTeg t G1 3 • 2  » 

•  •  UGE£,*2/N‘*.Tlt2.Gl3.2,‘  3AJT  »ILSS‘/ 

•  rT.‘Aa‘.TlS,‘ALTI“£TCR  SCM.C  FACT-R‘,'A*. 

•  ‘'A*(K»  ccns'anm.’^.su.a  , . •/ 

•  T7, T13, ‘STATIC  »»CSSJAC  C  OEFFCIE  NT*  .*  «*  • 

•  *:aso:«  ccssTALT‘,Tsr.Gi2.*,‘  fv«ft/«:ec,,2i‘/ 

•  t 7,‘sc‘,tis.*alt!»eter  l*g‘,t»«, 

•  ‘'Afio:«  c:nsta«.t‘,T62,;i2.»,*sccs‘7 


subroutine  tfaj 


74/74 


PTzl 


FTH  4.?»5S4 


11/00/42  15. 


1 


•J 


11 


15 


21 


25 


11 


35 


4] 


45 


*1 


55 


surf  cut:  vc  tsaji;»un,t.vf,vs,v«tj,«f,«s,*t®aj> 

c 

C  T4  AJ  C’CATES  AN  EXTERNAL  FLIGHT  >»3FILE  CF  OU®A*ICN 
C  641  SEC0V1S.IT  S  T»  4  TS  WT  T-4  *  LEViL  FLISHT  WITH  A  SPEED 

C  OF  401NILES  PE 4  Hi UP  FIR  5C1  5EC1V1S.it  CO-REvSES  "HE 

C  01  VE  FIR  10  SECS  WITH  DOWNWARD  VELOCITY  OF  GOCCF*/“IN, 

c  and  levels  if  at  hocft.this  trajectory  WAS  created  r0 

C  EVABLE  A  VEHICLE  TO  PERFO’h  A  TE«Clw  u®oate. 

c 

COHNCS  /HATH/OT»PI .HALFPI ,p» .twcpi.rpo 
CC  HHON/E  ASTH/CNEGA  .REQ.ES0.5EE 
C3«hCV/h6Tm/hT 

01 HESSION  AFlf.FI.  XS(NS).  xtrajinatji 

c 

C  SES “ENT  ONE 

IF<T-T1121. 22,22 

21  V2=0 . 

C  INI  UAL  HEISrtT  infeet 
HT  =1  10  00  . 

C  The  FOLLOWIvS  EQUATIONS  DESCM5E  A  OREAT  CIRCLE  »A'H 
C 

VHOR  =  SAHOOT*(REQ»H’l 
GAI4  r  VMC».«T-TO  l/IREO»Hr» 

CTLAT  S  SQST(CCS<GAH).C0'(SAh»»UC0S«tINCL)*SIN(GA»»)..2»» 
S  TL  »T  :  SIMTINCL)»S:V(SA*<> 

SALPhA  s  COS« f INCL»/CTlAT 
IF«GA"-GU)24,25.25 

25  IF1GAH-GI?- 26.26.24 

26  CAUPHASSQRTU-SALPHA.SALPHAI 
GO  T 1  27 

24  CALPha:-SQ"TU-SALPhA.SAL»hA) 

GO  TO  27 
C  SEG “ENT  TWC 

22  IF1T-TH40.40  ,41 

40  GO  TO  21 

41  IFIT-ITI .DELTA1I42.42.43 

42  V2s-50..(l.-CCS<OWE64<T-'t l»> 

GO  *0  44 

43  :F«T-«T2-0ELTA1I44.44,45 

44  V2=-100. 

GO  TC  44 

45  IFIT-T2146.46.24 

46  V2s-5D.4a.-C0S10NE«4<T-T211» 

GO  TO  44 

48  MT  *  UCC0.-130.4«T-Tll 

VhOR  =  GAHOCT* 1REB«HT1 
GA»  *  vhCP » 1 T-TO 1/ I REO»hT 1 

CTLAT  =  SQRT1C0S  1CAN1  .COS  (5A«I  .MC0S1T  IvCL>*I  INCGAWI  l.»2  1  I 
STL»T  s  S;MTIVCL>4=CV<51h1 
1  alpha  =  COSlTINCLl/CTLAT 
:F(GA«-61II3C. 31.31 

31  IF (GAH-G12  >32 . 32. 30 

32  CALPHA  =  51RTU  -IALPhA.S  AL»4A  » 

GO  th  27 

30  CAlphAi-SORTII-JAlPhA.SAlPhAI 

C  SEGMENT  THREE 
24  V2= 0. 
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subroutine  tpaj  T*n «  optpi  ftn  a.p.sg*  u/os/32  15.33 


Hi  i  1031. 

vhop  =  G**,oor««P£a»HT » 

S3  GAM  =  vhOc*(T-TO I/ISC0»«T> 

CTLAT  p  S0PT<C3S<GAM|.C;3<GA«l*<<C3S<TrNCL»*S:N<64M7».»2>» 

stlai  =  s:.n<t:ncl>«sin<gami 
SalPha  =  COS«T!NCL*/CIL*I 
:F(GA«-611133. 34,34 


41  34  :F«GA“-G12> 31.35.33 

35  CALPHAPSa«TU-SALPHA»$ALPHAl 

GO  rO  2 7 

33  CALPHA  =  -5aPTU-3*LPM*»34L9H*> 

C 

13  c  . ....oefini  rtoss  of  T»Aji;t:4r  va»iables . 

c 

C  «T4AJC1»  LATITUDE 

C  PT3Ad<2»  *  VELTCITT 

c  «t«aj«3»  t  »el::itt 

Tl  C  »T-AJt41  l  VE'.OCtTV 

C  XT  3  Adi  5  3  I  SPTCIMC  F3 3  CE 

c  *r3Aj(s>  r  specific  fo’ce 

C  XT3AdiT>  2  SPECIFIC  FO:CE 

C  xTPAd(57  UANOC»  ANGLE,  alpha 

53  C  XT3Ad<»*  HEIGHT 

C 

C  WHERE.  FGA  alpha  —  -30  DEG»EES»  X-f-Z  POINT  E-N-u  »ES  PE  C*  I  »ElV  . 
C 

27  VNJVhoR.CALPhA 

-5  VE  *  VHJS.salPha 

FE  =  -2. A  mOP.mJE.STLAI.CAlPMA 

FN  =  2«WHd».«IE.3TLAT.SALPHA 

FZ  :  GEE-IVHCP.WHCP/IAEO.HTn-Z.yHCP.yiC.CCSITINCLJ 

«t»Ad<n  a  aco:<ctla'» 

*?  »T»Ad«2»  P  VE 

X  TP  A  d  ( 31  P  VN 
(TPAd(41  —  VZ 
«r»Adl3»  p  FE 


«7«Ad<6»  —  FN 
-5  XT»Ad<7»  5  FZ 

«T»Ad<9»  P  HT 

»etu»n 

ENTPT  T»AdO 

C  HUIZONTAL  VEL.'ClTy  40  CPI LE S/H3U* 
lo:  V  H  OP  —  1013.37 

C  INI'IAL  latitude  IN  PAOIANS 
7  LA  r  p  0. 

C  TIME  OF  SEGMENTS  in  seconds 
TC  p  3. 

101  71  I  500. 

72  p  400. 

C  G3S  AT  CIPCLE  PA'H  I  N  CL  I N  AT  I  3  N  IN  >A3tANS 
TINCL  3  .7553345 
C  INITIAL  HEIGHT  IN  FEET 
111  M-  s  11037. 

GA«OCT  p  V"E4/f PEa^n'l 
Gil  —  -1.5TCT3T 
G 12  p  1.5  70  73  7 
G22  =  4.71233! 
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subroutine  r?»j 


7*  /  7  %  OPTsl 


UIC  S  .72^211^1*7 E-« 
OELr»=,31 
:«eo=Pi/OEL'» 

VAN  OCR  ANGLE  IN  SAOIANS 
*T9»J<  S»  *  -HALFP! 
P£7U«N 
ENO 
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sua* jurist  uspin 


7 «/?4  OPT  si 


f T N  4.F.5S* 


ll/0*/32  IS. 33 


11 


l  subroutine  uspin 

c 

C  USE  P-uP I  T  TE  N  SUBROUTINE 

c  this  routine  specifies  constats  »no  reaos 
5  c  :s  control  ano  statistical  o»r» 

c 

CC*«CN  /CCP2S/OALTS.DGRAVE.069AVN,DGA 
COH.HCN  /EAP-H/C-EGA.KEa.EIQ.GEE 
COHNCN  /CN*JL/Lc0a*,LINTC0,<«SS3 
10  CO«»ON  /«ATM/Q'TPt,MALF»! .PI .TWO®! ,P»0 

COHHt'N  /nOISOS/GTNOS(3I.ACnDS<3> 

COHhIN  /SI  GO /OTNANIC(10>.S3( 3 > , 3 - < S > .GFF < 3  I , G SF ( 3 > .G" < S >  . 

•  AB( 3 > , ASF ( ? > , A«( S » , 3A TO ,GDE .GO .6A , 

•  BASF.SPC.AL 

15  COH-CN  /StG“AS/8ASIGS.EDS!GS.0NSt6S.GA$IGS 

COF“IN  /VOA»o/C<l  >C*2tC«3 

c 

LOGICAL  LF08A.LINTOO 

c 

23  C 

NApELIST/CCN’PL/LFDBiC.LInTCO.RTSSO 

SAHEL  I  ST /S!GCS/0TNA**IC.68,6p  . G r F , GSF , G » . A 8 . A SF , A» , 

•  8  AP  0  *  G0EiGS0.GA.3ASF.jPC.AL 
SAsELlST/STArs/OALTS.DG»AVE.OG»A»N,OGA.GrSOS. ACN  OS  .  8  AS  I  G  S  .E  OS  I  GS  . 

25  •  oss:gs.gas:g3 

c 

c  SET  EARTi-AELATEO  CONSTANTS 
“EGA  *  .72«2U5l«7C*4 

-EQ  *  2.0P25GAE7 

53  C  ESI  =  1.0366PAJ177TS 

Esa  s  :.o 

GEE  s  12.2 

C 

c  set  vertical  l::p  oa»p:ng  constants 
35  C«l  =  3.3E-2 

CP 2  *  3.:3flTE-« 

CKJ  =  l.OE-S 
C 

C  COAPUTE  hath  ccsstants 
«3  PI  p  ACOSl-l.l 

ar’PI  =  0.2i*Pt 
halfpi  s  o.s;*pi 
twopi  s  2.c;«p: 

pp 0  =  PI/133. 

A3  C 

C  »EAO  ASO  ECiO  statistical  oata 

PEAO  (5.SIGCSI 
®EAO  <3. STATS) 

WRITEIGtSIGOI ) 

51  JRITE(i.lOO) 

UPtTE«6.STATS) 

WAITEI6.300) 

C 

r  9£A0  AND  ECHO  INITIALIZATION  CONTROL  PAJA-fEAS. 

55 

NOTE  THE  FOLLOWING  A  BCUT  CO'.T«OL  PATA<tep  k»SGO: 

C  <  XS  6  1  s  1  T  r  U  T  A  STATES  IMTIA.IZEO  IN  *SOCT?  8T  CVE®W»I'ING 
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SUBROUTINE  USMS 


74/7#  O^Trl 


U/!!/1?  15.33 


s 

. 

i 

.*  - 

i 

i 

*  »* 


63 


65 


73 


75 


-5 


=  ) 


C 

C 

C 

c 

c 


c 

c 

c 

c 


c 


xSO  I  NPUT  «ITh  ONE-SISNA  VA.UCS  F*0«  NA-ELIST  *StG3'*. 
FxSGO  =  2  T’UM  STATES  INt'IA.IZEO  87  "ONTE  CA»LO  P*CCESS 
USING  ONE-SIG'A  VALJES  F=0*  NA*ELl$T  "SIGeS*. 

otherwise  all  r  =  urn  states  initializes  facn 

XSO  INPUT  « 5T  1-t  ir;«  15  TABLE  *-2  "F  SCFE  *ANUAL). 

3EAQ  <S.C0N'-L1 
V8  I T  E<  6 «  C  ON"  =  L  > 
w*!TEI6.«C3> 

PAINT  TABLES  OF  STATISTICAL  I  N  P  U  T 
CALL  STABLE 

CON  V£P  T  INPUT  JATA  TO  CCHPUTATI oval  UN  I  * S  AND  ECHO  IT  ONCE  *0*E 
CALL  NUUNIT 
y®!TE«S.SIGC:» 

WAI T  E( 6t  S  TATS) 


P  E  T  Uc  N 

i c 3  f:phat«*c*,5«.*note:  the  sigcs  namelist  con*ains  sig*a  • 

•  •values  for  initialization  :f  tajth  «ooel  states  c»s:i  unoe»  • 

•  •cos‘®ol  of  kxeoo  as  foll:ws:*//tij,*kxsco"I7,*ac*ion» 

•  /MS,*l*.T2J,*x30  SET  TO  0N£*SIG»A  VALUES  F»0«  SIGOS*. 

•  /TI5,*2*,T2J,*IS3  «CNTE  CARL0E3.  EACH  ''F  NS  SAHPLES  IS  GAUSSIAN* 
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